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ABSTRACT 

We investigate statistically the dynamical consequences of cosmological fluxes of 
matter and related moments on progenitors of today's dark matter haloes. These 
haloes are described as open collisionless systems which do not undergo strong in¬ 
teractions anymore. Their dynamics is described via canonical perturbation theory 
which accounts for two types of perturbations: the tidal field corresponding to fly¬ 
bys and accretion of dark matter through the halo's outer boundary. 

The non-linear evolution of both the entering flux and the particles of the halo 
is followed perturbatively. The dynamical equations are solved linearly, order by or¬ 
der, projecting on a biorthogonal basis to consistently satisfy the field equation. Since 
our perturbative solution of the Boltzmarm Poisson is explicit, we obtain, as a result, 
expressions for the N-point correlation function of the halo's response to the pertur¬ 
bative environment. It allows statistical predictions for the ensemble distribution of 
the inner dynamical features of haloes. We demonstrate the feasibility of the imple¬ 
mentation via a simple example in the appendix. We argue that the fluid description 
accounts for the dynamical drag and the tidal stripping of incoming structures. We 
discuss the realm of non-linear problems which could be addressed statistically by 
such a theory, such as differential dynamical friction, tidal stripping and the self grav¬ 
ity of objects within the virial sphere. 

The secular evolution of open galactic haloes is investigated: we derive the ki¬ 
netic equation which governs the quasi-linear evolution of dark matter profile in¬ 
duced by infall and its corresponding gravitational correlations. This yields a Fokker 
Planck-like equation for the angle-averaged underlying distribution function. This 
equation has an explicit source term accounting for the net infall through the virial 
sphere. Under the assumption of ergodicity we then relate the corresponding source, 
drift and diffusion coefficients for the ensemble-average distribution to the underly¬ 
ing cosmic two-point statistics of the infall and discuss possible applications. 

The internal d5mamics of sub-structures within galactic haloes (distortion, 
clumps as traced by Xray emissivity, weak lensing, dark matter annihilation, tidal 
streams ..), and the implication for the disk (spiral structure, warp etc... ) are then 
briefly discussed. We show how this theory could be used to (i) observationally con¬ 
strain the statistical nature of the infall (ii) predict the observed distribution and cor¬ 
relations of substructures in upcoming surveys, (iii) predict the past evolution of the 
observed distribution of clumps, and finally (iv) weight the relative importance of 
the intrinsic (via the unperturbed distribution function) and external (tidal and/ or 
infall) influence of the environment in determining the fate of galaxies. We stress that 
our theory describes the perturbed distribution function (mean profile removed) di¬ 
rectly in phase space. 

Key words: Galaxies: haloes, kinematics and dynamics, statistics, Cosmology: Dark 
Matter 
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1 INTRODUCTION 

It now appears clearly that the dynamical (azimuthal in¬ 
stabilities, warps, accretion), physical (heating, cooling) and 
chemical (metal-poor cold gas fluxes) evolution of galaxies 
are processes which are partly driven by the boundary condi¬ 
tions imposed by their cosmological environment. It is there¬ 
fore of prime importance to formulate the effects of such an 
interaction in a unified framework. 

Modern digital all-sky surveys, such as the SDSS, 
2MASS or the 2dF provide for the first time the opportu¬ 
nity to build statistically relevant constraints on the dynam¬ 
ical states of galaxies which can be used as observational in¬ 
put. Other projects, like Gaia or Planck, will provide small- 
scale information on our Galaxy and its environment and 
will soon allow detailed confrontation of the predictions of 
models with the observations. We ought to be able to draw 
conclusions on the internal dynamics of the halo and its inner 
components and constrain their statistical properties. 

Unfortunately, it is difficult to study the response of 
haloes to moderate amplitude perturbations. Current N- 
body techniques suffer from resolution limitations (due 
to particle numbe r and drift in orbit integration, see e.g. 
IPnwer et alJ J2Q^, iRinnevI i2nn4h . for a discussion of such 
effects) that hide to some extent linea r collective effect s 
which domin ate the response of the halo iWeinbergI il998bh . 
|m uralil h999l) ')^. Simulations on galactic scales are also of¬ 
ten carried without any attempt to represent the cosmo¬ 
logical variety arising from the possible boundary condi¬ 
tions (the so-called cosmic variance problem). This is be¬ 
cause the dynamical range required to describe both the en¬ 
vironment and the inner structure is considerable, and can 
only be ach i eved f or a limited number of simulations (e. g. 
iKnebe et"^ i2004h . iGill et all i2004h . IPiemand et alJ j2004h '). 
By contrast, the method presented below circumvents this 
difficulty while relying on an explicit treatment of the inner 
dynamics of the halo, in the perturbative regime. Specifically, 
our purpose is fo develop a tool to study the dynamics of an 
open sfellar system and apply it to the dynamic of a halo 
which is embedded into its cosmological environment. One 
can think of this project as an attempt to produce a semi- 
analytic explicit re-simulation tool, in the spirit of what is 
done in N-body simulations with zoomed-in initial condi¬ 
tions. 

The concept of an initial power spectrum describing 
the statistical properties of the gravitational pertu rbations 
has p r oved very useful in cosm ological studies (e.g. lPeeblesI 
il98(1i . lRernardeau etalJ i2n02h 'l. The imderlying paradigm, 
that gravity drives cosmic evolution, is likely to be a good 
description at the megaparsec scale. We show below that a 
similar approach to galactic haloes is still acceptable, and 
marginally within the reach of our modeling capabilities. The 
description of the boundary is significantly more complex, 
but the inner dynamics of hof components is better behaved. 
Here, we describe a stable system which imdergoes small in- 


^ it has been argued that shadowing jEarn &: Tremainel Jl99]h ') will 
in practice allow for another orbit to correct for the drift, but this is of 
no help to resonant processes because it requires that the same orbit 
does not diffuse for a few libration periods. 


teractions, rather than an unstable system in comoving coor¬ 
dinates undergoing catastrophic multi-scale collapse. 

The purpose of this investigation is to derive analytically 
the dynamical response of a galactic halo, induced by its (rel¬ 
atively weak) interaction with its near environment. Inter¬ 
action should be understood in a general sense and involve 
tidal potential interactions (like that corresponding to a satel¬ 
lite orbiting around the galaxy), or an infall where an external 
quantity (virialized or not) is advected into the galactic halo. 

With a suitable formalism, we derive fhe propagation 
of an external perturbation from the near galactic environ¬ 
ment down to the scale of fhe galactic disk through the 
dark matter halo. We essentially solve the coupled collision¬ 
less Boltzmann-Poisson equations as a Dirichlet initial value 
problem to determine the response of fhe halo to infall and 
tidal field. The basis over which the response is projected can 
be customized to, say, the universal profile of dark matter 
haloes, which makes it possible to consistently and efficiently 
solve the coupled dynamical and field equations , so long as 
the entering fluxes of dark matter amounts to a small pertur¬ 
bation in mass compared to the underlying equilibrium. 

In a pair of companion papers, lAubert fe PichorJ 
<200.5allJ) described the statistical properties of fhe infalling 
distribution of dark matter at the virial radius, R 200 a 
function of cosmic time between redshift z = 1 and today. 
These papers focused on a description of fhe one- and two- 
point statistics of fhe infall towards well formed L* dark mat¬ 
ter haloes. All measurements were carried for 15 000 haloes 
undergoing minor mergers. The two-point correlations were 
measured both angularly and temporally for the flux densi¬ 
ties, and over the whole 5D phase space for the expansion 
coefficients of fhe source. 

Together with the measurements presented there, we 
show in this paper that the formalism described below 
will allow astronomers to address globally and coherently 
dynamical issues on galactic scales. Most importantly it 
will allow them to tackle problems in a statistically rep¬ 
resentative manner. This investigation has a broad field 
of possible applications. Galaxies are subject to bound¬ 
ary conditions that reflect motions on larger scales and 
their dynamics may constrain the cosmology through the 
rate of merging events for example, or the mass distri¬ 
bution of satellites. Halo transmission and amplification 
also fosters com municat i on bet ween spatially separated re¬ 
gions, (see e.g. iMuralil il999h ') and continuously excites 
the disk structure. For example, spirals can be induced 
by encounters with satellites and/or by mass i niection 
te.g. lToomre & ToomreHl973l . lHoward & BvrdHl99ril) '). while 
warps resul ts from torque interactions wi t h the surround¬ 
ing m atter ihopez-Corredoira et alJ i2002h . Iliang & BinnevI 
jl2^). Therefore the proportion of spirals and warps con¬ 
tains information on the structure's formation and environ¬ 
ment. The statistical link between the inner properties of 
galactic haloes, and their cosmic boundary can be reversed 
to attempt and constrain the nature of the infall while inves¬ 
tigating the one and two point statistics of fhe induced per¬ 
turbations. This is best done by transposing down to galac¬ 
tic scales the classical cosmic probes for the large-scale struc¬ 
tures (lensing, SZ, etc... ) which have been used successfully 
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Figure 1. The action angle, (I, w), - spherical coordinate, (r, v), trans¬ 
formation. The dark matter particle at running spherical coordinate 
{r,d,(p) describes a rosette in the orbital plane orthogonal to its mo¬ 
mentum, L. The line of node of the orbital plane intersects the x — y 
plane at a constant (in spherical symmetry) angle with respect to 
the X axis. The orbital plane is at an angle fi = acos(Lz/L) to the 
X — y plane. The particle polar coordinate in the orbital plane w.r.t. 
this line of node is ip. The angle coordinates, W 2 , is measured along 
Ip but varies linearly with time by construction. Finally the radial an¬ 
gle ivi varies with rad ius between peri-apse and ap o-apse. (Strongly 
inspired from Fig 1 of iTremaine & Weinberelli984ll '). 

to characterize the power-spectrum of fluctuations on larger 
scales. 

The outline of this paper is the following: we describe 
in Section 13 the linear response of a spherical halo which 
undergoes cosmological infall of dark matter, and compute 
the induced correlations in the inner halo; Section|3presents 
the second-order perturbative response of the galactic halo to 
the infalling flux; ( Appendixlolgives the higher order correc¬ 
tions to the dynamics and addresses the issue of dynamical 
friction). Section 0] derives the Fokker-Planck equation that 
the cosmic mean halo profile obeys in such an open environ¬ 
ment. Section|3describes briefly possible astrophysical appli¬ 
cations. In particular, it is discussed how the statistical anal¬ 
ysis of mean and variance properties of galactic haloes and 
galaxies can be compared to the quantitative prediction of 
the concordant ACDM cosmogony on those scales. We also 
show how to revert in time observed tidal features within our 
Galaxy, or in external galaxies. The last section draws conclu¬ 
sions and discusses prospects for future work. 


2 THE SPHERICAL HALO: LINEAR RESPONSE 


IVauterin & Deionghel Jl99(j) . iBertin et alJ il994l by adding a 
source term to the collisionless Boltzmann equation. Since 
the formalism is ofherwise fairly standard, we will present it 
relatively swiftly. In a nutshell, the dynamical equations are 
solved linearly while projecting over a biorth ogonal basis to 
consistently satisfy the Pois son equation (e.g. iKalnai j h97lh 
iKalnaifi il97(jll^lnaisl il977h ). The dynamical equation of an 
open system characterized by its distribution function, F, to¬ 
gether with the field equation, read formally: 

afF+jH-l-tp^F} = 5'", and = TttS ^ dvF, (1) 

where H is the Hamiltonian of the system, { } is the usual 
Poisson bracket and stands for the potential, the 

perturbing exterior potential and incoming source term. The 
latter, s‘^(r,v, f) accounts for the entering dark matter at the 
viria l radius and is discu ssed in detail below (Sectionl2.3l see 
alsn lAiibert et alJ i2nfl4ft and lAubert fe PichnrJ In a 

somewhat unconyentional manner, t) refers here to the 
external potential, i.e. the tidal potential created by the per¬ 
turbations outside the outer boundary of the halo (, i.e. R20o)- 
Let us expand the Hamiltonian and the distribution 
function, F, as: 

¥ = F + £f, and H = Hq + eip, (2) 

where we assume that eyerywhere in phase space 
£ = m/M 1 i.e. that the mass of the perturbation, m, 
is small compared to the mass, M, of fhe unperturbed 
halo. In Eq. Q, / represents the small response to the 
perturbations, F represents the equilibrium state and ip the 
small response in potential. Putting Eq. j3 into Eq. ^ and 
reordering in e yields the linearized Boltzmann equation : 

dt dl dw ^ dw dw ^ dl 

where I and w are conjugate canonical yariables which are 
described in the following section. 

2.1 The Boltzmann equation in action-angle 

The most adequate representation of multiply-periodic inte- 
grable systems relies on the action-angle yariables^, since res¬ 
onant processes will dominate the response of the liye halo, 
and are best expressed in those yariables. We will use yec- 
tor notation for simplicity. The details of fhe computation of 
th ese yariables i s discussed in Appendix IbI following work 
bylM ura ilHli (see also Fig. 0). This achieyes separation 
of yariable between the phase space canonical yariables (an¬ 
gle and actions) on the one hand, and time on the other hand. 
We denote as usual the set of action yariables by I and angle 
yariables by w (see Appendix|Bj. The rates of change of an¬ 
gles is cu = dw/dt. Along the multi-periodic orbits, any field, 
Z, can be Fourier-expanded with respect to the angles as: 

Z(r,v,t) = £Zk(I,t)exp(!k-w). (4) 

k 


In the following section, we extend to open spheri¬ 
cal st ellar systems the fo rmalism d e yelop ed by iKalnaisI 
il97(jl (for stellar disks), lAo HeL^ ll9Z^ (for gaseous 


disksl . lFridman fe Poliacbenknl It 98^ TVgimin^^^^inbgrgl 
h984ll and e.g. IPaimeT&Papak^oulTwS^riMur^ h999ft . 


Conyersely 

^ Note that (w, I) are canonical variables, and as such preclude noth¬ 
ing about the evolution of the system. They simplify the expression 
of the linearized equations, order by order. 
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^/dwexp( 


l\^ ■ w)Z(r, V, t), 


(5) 


where k = {ki,k 2 ,k^) is the Fourier triple index conjugate to 
the three angles w = (zv-i, W 2 , W 3 ) . Given Eqs. Q-O, the lin¬ 
earized Boltzmann equation in such a representation is: 

^ +*k' cu/kjl, t) = zk- — t) + f[(l, t)) +sl(I,t).( 6 ) 


Here tp is the potential perturbation created by the halo's 
inner density fluctuations and ip‘^ the potential perturbation 
created by external fly byes. The gravitational field of incom¬ 
ing particles is accounted for by the source term s'^. The solu¬ 
tion to Eq. may then be written as: 


fkilt) = [ exp{ik-w{T-t))x 

J —00 


dF 


'k • [ipkil t) + fiil t)] -f t) 


dr. 


(7) 


Eq. 0 assumes that the perturbation has been switched on a 
long time ago in the past so that all transients have damped 
out.^ 


2.2 Self-consistency 


Eq.O can be integrated over velocities and summed over k 
to get the density perturbation: 


F(r- 


A) = L dr / dv exp(ik • a7(T — f) + !k ■ w) 

J —CO J 


dF 


df ['/’k(kT) +'/’k(kT)] +4(kT) 


( 8 ) 


Let us expand the potential and the density over a bi- 
orthogonal complete set of basis functions such that: 

= p(r,t) = X]fln(f)pt"l(r), (9) 

n n 


= 47rGpN , J !/;H*(r)p[Pl(r)dr = J", 


( 10 ) 


(where is the complex conjugate of !/^["l(r)). We nafu- 

rally e xpand fhe exfernal potential on the same basis iKalnaisI 
<197lh l as: 


n 


(11) 


Thus, the coefficients fln are representative of fhe density and 
potential perturbations in the halo itself, at r < F 200 / while 
the coefficients fen represent the potential created in the halo 
by density fluctuations at r > F 2 oo- Taking advantage of bi¬ 
orthogonality Eq. 13 is multiplied by tpp(r) and integrated 
over r, which yields: 


(0 = zf dr ff dvdrexpjzk ■ (ujT — t) + !k ■ w)i/’[Pl*(r) X 
1, J —00 JJ 


k 

E*k- ^ [«n(T) +fen(T)]l/>[‘’'(I) -FS^(I,t) 


( 12 ) 


^ Mathematically, we only retain the particular solution to Eq. |3, 
while assuming that the homogeneous solution did not hit long- 
lived resonances. 


We may now swap from position-velocity to action-angle 
variables. Since this transformation is canonical dvdr = 
dwdl. In Eq. tl2l only ( r) depends on w, so we may carry 

the w integration over which yields !/)|^^*(I). Eq. 
then becomes : 

flp(f) = (27 t)^E f dr f dlexp(!k • wjT — t)) X (13) 

E'k- ^ [fln(T) +fen(T)]tp[P’*(I)l/4"'(I) +Sk(kT)l/'Jf'*(I) . 

At this point, it seems natural to expand the source term on a 
basis too, but unlike the previous one, this basis should also 
describe velocity space. We admit for now fhaf such a basis 
(p„{i,v) exisfs, and write: 

s''(r,v,t) = Ecn(f)i^>t"'(rw), so (14) 

n 

Sk(kT) = (15) 

n 

where i7^"^‘^(I) is the angle transform of (r, v) (see Eq. i24l 
below). The coefficients Cn are representative of the mass 
exchange between the halo and the external world. The 
sum in Eq. iTst spans velocity space as well as configura¬ 
tion space, and therefore involve significantiy more terms . 
Such an expansion is performed in I Aubert fe PichonI i20n,3ah 
to constrain the source function measured in cosmological 
simulations. Calling a(T) = [apr),-■ ■ ,a-aA) ■ ■ ■], = 

[fei(T),-• •,fen(T) ■ • ■], and c(t) = [ci(t), • ■ ■ , Cn(T) ■ • ■], we 
define two matrices, K and Q. The matrix K has elements 
f^p,n defined by: 

J<p,n(T) = (27T)3WdIexp(zk-wT)zk-^!/;[Pi*(I)!/:j^‘’'(I),(16) 

which depend only on the halo equilibrium state. The matrix 
Q has elements 

Qp,n(T) = {2nfYl /dlexp(!k • wt) 4"'"(I) i/'lf'*(1), (17) 

k J 

which depend only on the source's expansion basis. Equation 
iTsl then becomes: 

a(t) = f dT (K(t - t) ■ [a(T) -f b(T)] -f Q(t - f) • c(t)) . (18) 
J — 00 

The kernels K and Q are functions of the equilibrium state 
distribution function, F, and of the two bases, cj)l'^l(r, v), and 
!p["l(r) only. They may be computed once and for all for a 
given equilibrium model. Assuming linearity and knowing 
K and Q, one can see that the properties of the environ¬ 
ments (represented by b and c) are propagated to the inner 
dynamical properties of collisionless systems (described by 
a). We may perform a "half" Fourier transform with respect 
to time. In the limit where the transients may be neglected, 
which implies that the system should be stable, this trans¬ 
form amounts to a Laplace transform with p = icv + . 

Temporal convolutions are then replaced by matrix multipli¬ 
cations and Eq. <181 becomes: 

a(a;) = (1 - K(a;))^i • [K(a;) • b(a;) + Q(a;) • c(a;)j . (19) 
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Figure 2. A typical rosette orbit in its orbital plane; the intersection 
with the R 200 sphere is shown, together with the corresponding ve¬ 
locity vectors, both entering and exiting. The net flux of such quan¬ 
tities enters Eq. 12.51 and characterize the source of infall perturbing 
the halo. Note that by construction, in the linear regime all infalling 
material re-exits R 20 O/ since the perturbation evolves along the un¬ 
perturbed orbits. This is to be contrasted to the situation presented 
in Fig. 13 where dynamical friction is qualitatively accounted for. 


and a similar expression for Q[ 2 ],n(r/T), involving the 
source's expansion basis, the mass flux may be written as a 
convolution: 

|Ov(r,f) = /dT(K 2 (r,T-f) ■ [a(T) + b(T)] + Q 2 (r, t - f) • c(t)) . 

J —CO 

After half Fourier transforming wifh respect to time, we get 
pv = K[ 2 ] (r, w) • (1 - • [b + Q • c] + Qpj (r, w) • c. (22) 

We will return to Eq. I22> in Section|3 


2.3 Sinks, sources and tidal field 

Let us now turn to an explicit description of fhe source term 
(s*^, hence c), and the tidal field hence b) entering Eq. 

We consider here a source at the virial radius correspond¬ 
ing to cosmic infall. Note however that we might have con¬ 
sidered just as well sinks reflecting the presence of a super 
massive black hole at the center of the host galaxy or the de¬ 
flection/absorption of orbits due to a galactic disk. 


2.3.1 Source of infall at R 200 

A possible Ansatz for fhe source term consistent with the first 
two velo city moments of the entering matter has been pro¬ 
posed bv lAubert et alJ i2004h . Following them s‘^(r, v, t) can 
be written as 


In this expression, 1 is the identity matrix, and K and Q in¬ 
clude Heaviside functi ons before Fourier transform to ac¬ 
count for causality fsee lAubert et alJ i2004h for details). Sec- 
tion lDSI gives an explicit expression for K(a;). Nofe fhe differ¬ 
ence between w, the angular frequency of the orbits, defined 
above Eq. and iv, fhe half Eourier fransform variable as¬ 
sociated with time which appears in Eq. <191 . Here b and c 
could be given deterministic functions of t ime, or stochastic 
random fields (characterized statistically in lAuhert fe Pichn^ 
<200.5ah l. In contrast, a describes the detailed response of fhe 
halo in phase space within 1^200 ■ 


2.2.1 Higher moments 

The second moment is obtained by multiplying j7| by v and 
by performing an integration over velocities. Summing over 
k leads to: 


|Ov(r, 


', t) = dr dvexp(!k • a;(T — t)-f ik • w) 

V J —00 J 


dF 


'k ■ -^v t) + flih t)] -f vs^(I, t) 


( 20 ) 


Using the same biorthonormal expansion as above, we may 
express the mass flux as a function of fhe coefficients On and 
bn (associated with the potential perturbations of external 
origin). If we define the following new tensors: 

%,n(TT) = X] [ dvexp(!k-wT-f!k-w)v !k- ^!/;["1(I),(21) 
k 


Cn(f)S'ii:(’'')'''m(El) J , 

where m stands for the two harmonic numbers, {£,in) and 
Ymjn) = y™(n) is the usual spherical harmonic. The Dirac 
function — R 200 ) appears because the source terms are 
located at the virial radius in our representation'^. This equa¬ 
tion corresponds to the parameterization of as: 

<fM(r,v) = ya(v)ym(n)bD(r-R 20 o). 

= g«(t’)ym(n)ym'(r)bD(?--k200), (23) 

of Gaussian functions, ga, covering the radial velocity com¬ 
ponent and spherical harmonics for the angle distribution, T, 
of the v elocity vector and orient ation, n = { 6 ,ip) of the in¬ 
fall fsee lAubert fe PichorJ 1200, 5ah for defails). Here we have 
n = [m, a] = [m, m', a] = \£, m, t', m', a]. Erom Eq. <1,5< 

j d^wexp(-zk-w)(^["'(r,v). (24) 

Wifh Ea.l23l. Eq. <241 becomes 

= ■^^^y'd^wexp(-zk-w)y„[n(I,w)] X 

ya(v[I,w])<5D(r(I,w) - ^ 200 ) • (25) 

^ This choice is ma inly j ustified by the measurme nts performed in 
lAubert et alJ 1200411 and lAubert &: Piehorj l2005ah and stands as a 
good compromise between a relaxed halo inside this boundary and 
a low contribution of the orbits of relaxed particles to the flux. 


s''ii,v,t)=Y^Yra{Cl)Suir-R2oo) I] 
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We can make use of the function occurring in Eq. i25\ since 
Wy = Wr(r, I) (given by Eq. Therefore Eq. <251 reads: 


J dK.exp{-;kw)ym[0{I,w)l 




d^w 
(2 

ga(v[I,w]) 
d^w 


-—SuiWr-Wr[R200,i\), (26) 


(277) 


\dWr/dr\^^ 

3 exp(-!k'W)ym[n(I,w,a7p[R200/I])] x 


ga(v[I,W,rop(R 200 /I)]) 


COrjl) 

|f(^200xl)l 


exp(-ifcr-H;r[R 200 xI]) • 


In Eq. <26> we sum over all intersections of the orbit I with the 
R 200 sphere, at the radial phase corresponding to that inter¬ 
section with a weight corresponding to a;r/|f| (see Eig. j3). 
Note that Eq. <26> involves d^w = dzv 2 dzv 2 . 


2.3.2 Tidal excitation from beyond R 20 Q 

The tidal potential is given as a boundary condition on the 
virial sphere and deprojected in volume. Let us call 
the harmonic coefficients of the expansion of the external po¬ 
tential on the virial sphere. We expand the potential over the 
biorthogonal basis, (uff, d^'™) (see Appendix|Bj, so that 

r{r,n,t) = ^ b'^(t)y'"(n) 

n,e,m V -'’^200 / 

= £fon(t)V-["l(0. (27) 

n 

where (r) = y^’(n)M('"(r). The first equality in Eq. <27> 
corresponds to the inner solution of the three-dimensional 
potential whose boundary condition is given by Y^{Cl)b'^^ 
on the sphere of radius R 200 (defined below). Since the basis 
is biorthogonal, it follows that 



It is therefore straightforward to recover the coefficient of the 
3D external potential from that of the potential on the sphere. 


2.4 Induced correlations in the halo 

Our purpose is to characterize statistically the response of 
the dark matter halo to tidal perturbation and infall. This 
is besf done by computing the N-point statistics of fhe per¬ 
turbed density field. Let us start with the two-point corre- 
lahon. From Ea.iT^ the variance-covariance matrix of fhe 
response is given by 

(Ja-r^^ =([K.b + Q-c]-(l-K)-i. 

(l-K)^i*T • [K-b + Q-c]^*). (29) 

This expression of fhe n x n matrix, (a • a*^) involves 
autocorrelation terms like the components of (b ■ b*^) (the 
tidal field) and ( c • c*^) (the source of infall), but also cross¬ 
correlation terms such as the components of (b • e*^). For a 
spherical harmonic basis, the induced density perturbation 
reads (see Eq. <B2> in Appendix I bI 


p(r,n,t) =^flnp„(r) = a'}m{t)yr{(^)d'emP) , (30) 

n n^m 

The functions (r) depend on the chosen basis. An exam¬ 
ple is given by Eq. <B3> . Again, n stands here for n,£,m, 
respectively the radial and the two angular 'quantum num¬ 
bers'. As a consequence the two-point correlation function 
for the perturbed density reads 

(p(r, n + AO, t + At)p(r', Cl,t)) = ^ ^“(n) x 

nimn'i'm' 

Yf*(n + An)dUr)dt,yn,{r’){aUt)4m'it + Af)). (31) 

The statistical averages, {a"A^)a”',^,{t + At)) are given by the 
temporal inverse Fourier transform of Eq. <29t . If fhe per¬ 
turbation is stationary and statistically rotationally invariant 
(fl",„(f)fl”/’„,/(f + At)) =Q"'(At)<5”''(5|'. The correlation func¬ 
tion then obeys 

(p(r, n -f AO, t At)p{r', Ci,t)') = 

£ PdcosA))dUr)dt{r')Cr'm , (32) 

nn'£m 

where 7 stands for the angle between fl and Cl'. Evaluat¬ 
ing Eq. <321 for 7 = 0, At = 0, r = r' gives a measure of 
fhe cosmic variance of fhe amplitude of the response of the 
halo as a function of radius r. The full-width half maximum 
(FWHM) of (p(r, n + An, t)p{r, Cl, t)) is a measure as a func¬ 
tion of time, t, and radius, r, of fhe angular extent of the en¬ 
semble average mean polarization. Conversely, the FWHM 
of (p(r, Cl, t)p{r -f Ar, Cl, t)) is a measure of its radial extent 
in the direction Cl. Note that Eq. j7| together with <18t yield 
a description of the response both in position and velocity. 
For instance, Eq. allow us to predict the induced corre¬ 
lations amongst streams. Applications of Eqs. <29l-i3^ (and 
their non-linear generalization in Section I3.2> will be dis¬ 
cussed in greater details in Section The actual implemen¬ 
tation of Eqs. <29> -f??ls carried in a simplified framework in 
SectioniBll 


2.4.1 Link with propagators 

Lef us emphasize that the splitting of the gravitational field 
into two components, one originating outside of l? 200 x and 
one from the inside, via point particles obeying the distribu¬ 
tion s®(r, V, t) is somewhat ad hoc from the point of view of 
the linear dynamics. It is convenient from the point of view 
of fhe measurements, and crucial for the non-linear evolution 
(described below), or the ensemble average, as shown above 
It allows us to specify the statistical characteristics of the 
infall without having to refer to the properties of fhe object 
on which this infall occurs. 

We discuss in appendix ^ the formulation of fhe re¬ 
sponse of a self-gravitating sphere in terms of a propaga¬ 
tor (i.e. the Green function of the collisionless Boltzmann- 
Poisson equation). 


® One should account for the fact that tpy should be switched on 
long before any particles enter R 200 since no particle is created at 
the boundary. 
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This formulation is mathematically equivalent to the approach described above, but there we relied on Gauss's theorem to 
reproject all the information beyond R 200 back onto the virial sphere. This information involves two contributions: one relative 
to particles beyond R 20 O/ which contribute to the tidal field, the other relative to particles entering K 200 which contribute to 
5*^(0, v)) 

The main asset of this formulation is to localize the boundary, which is possible since the interaction is purely gravitational, 
at the expense of having two sources of different nature. In particular, this implies that the environment may be characterized 
once and for all, independently of the detailed nature of the inner halo. 

In this section, we assumed that the polarisation of the halo was linear. This hopefully provided some insight for some 
aspects of the dynamics, but effectively ignores non-linear phenomena such as dynamical friction or tidal stripping. Let us now 
expand perturbation theory to higher orders. 


3 NON-LINEAR PERTURBATIVE RESPONSE 

In the following, we describe, using perturbation theory, the non linear response of the halo to material entering at the virial 
sphere. It is assumed that the perturbation is first-order in the hierarchy, and that the halo is dynamically stable. This should 
warrant the validity of the expansion. We use the angle-action variables of the unperturbed system as canonical variables and 
investigate the non-linear evolution of the infall and the tidal excitation. 

In essence, the key is to expand the potential onto the biorthogonal potential density basis which allows us to decouple position- 
velocity and time (i.e. perform a separation of variable), and solve in turns each order of the perturbation expansion. 


3.1 Perturbative expansion 

Recall that the dynamical equation of an open system characterized by its distribution function, F is given by Eq. Q. Let us 
expand again F as 

F = F + and H = Hq- f (33) 

n n 


where the unperturbed equilibrium is characterized by the distribution function, F (I). Note that (n), the order of the expansion 
should not be confused with n = {n,£,m). Finally, it is assumed that the external perturbation enters as a first-order only, i.e. 
s‘^ oi e and xp'^ cx e. In short, the rewriting of Eq. 0 to order yields: 


dt 




n-1 


• 'k[ ^ + 4®^ 

k=l 


(34) 


In the following, we solve Eq. <341 recursively, order by order, to recover the perturbative response of the halo to the tidal 
interaction and infall. We expand both the potentials and the source term over a biorthogonal basis, so that, with referring to 
the order in the hierarchy and to the label in the basis 


P P 


4(kf) = Ecp(o4‘’'w ■ 

p 


(35) 


Recall also that the superscript, IPl, in Eq. <3.5> spans discretely a 3D or 5D space depending on the type of function basis. The 
first-order solution for flp was given Eq. <131 . Let us turn to the higher order equations. 


3.1.1 Second-order perturbation theory 

( 2 ) 

The second-order equation for flp ' reads 

af\t) = {InfY. J j • !kexp(!k • a7[T - t]) j -f 

( 27 T)^y drEy dlexp(!k- w[t- t]){/W(T,w,I),</-W(T,w,I)}_^,/;lfl*(I), 
where {f ^^\is the Poisson bracket of the perturbation to first-order. Now for a set (/, xp) we have 
{fH}k = j dwexp(-!k-w) |E/ki(I)exp(!ki •w),E!/’k2(I)exp(!k2'W)| . 


(36) 


(37) 


Therefore 
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ki+k2=k V / ki+k2=k 

where the sum is over kj with k 2 = k — kj. Given Eq. 0, Eq. <361 may be rearranged as 


( 38 ) 


apht) = J ‘^'^1 «q?('^i) |^E J dl^/'k''(I)'/'k'*(1) ' ;kexp(ik • w[ti - f]) j + 

(2nf [ dTi f dT 2 f dIexp(;k-a7[Ti-t]) ^ + fcqj( ti)] x 

J J V k ^ 


E 

ki+k2=k 


exp(!ki • a;[T 2 - ti]) 


• iki [42V2) + fcq2(T2)]!/’kJ"'(I) +Cq2[T2]‘7-^f‘*''(I) 




[qi] 


f't 


(39) 


Note that the r.h.s. of Eg. <391 is linear in while it is quadratic in involving products such as a a, a c, b c and so 

on. More generally, the perturbation theory at order (w) is linear in Note also that Eq. involves a double ordered time 
integral over Tj and T 2 of the source coefficient, Cqj (tj) and Cq^ (T 2 ), which accounts for the fact that, non-linearly, the relative 
phase of the accretion events matter (Eq. <0121 gives the analog to Eq. <391 in the complex frequency plane). Eq i3^ includes in 
particular a term like 




exp(!(ki + k2) • a;[Ti - t]) exp(!ki • w[t 2 - ^ —^'/’k?'-(“qlVi) + V ('^1)) % ('^2) 


(40) 


qi,q 2 


which involves the rate of change of the source term with respect to action variation (via / 81) modulated twice over time 

as exp(!(ki +k 2 ) • w[ti — f])exp(!ki • a7[T2 — Tj]). 

The second order solution can be synthetically written by introducing tensors K 2 and Q 2 similar to those defined in Eq. 
and Eq. ffTi to express the first order solution as Eq. <181 . These latter tensors will now be referred to as Kj and Qi. Specifically, 
the components of these tensors are defined as: 


8 f 


■!k. 


(Ki)p,qiKi -t] = (K)p,qi[Ti -t] = {Inf^ [ dlexp(!k- w[ti - t])i/)Jf'* 

k J 

(K 2 )p,qi,q 2 [Tl - hT 2 - Ti] = (27r)3 ^ f dlexp(zk • w[ti - f]) ^ exp(!ki • w[T 2 - Ti])|^ • iki 

c •’ ki+k2=k If 




(41) 

(42) 


while Qi involves replacing 8F/8I • k by ■ For instance. 


(Q 2 )p,qi,q 2 [Ti-hT2-Ti] = (27r)3X] / dlexp(zk • w[ti - t]) X] |exp(!ki • w[T 2 -Ti])n^'^'‘’^',i/;j^‘l^^'| 

V ki+k 2 =k 


This implies in particular that Qi = Q given by Eq. <171 . Note that each component of K 2 has the same complexity as Ki, i.e. the 
perturbation theory is linear order by order; on the other hand it involves all the couplings in configuration space, hence the 
double sum in k. With these definitions, Eqs. iTsl and <391 read formally 


= Ki • + b] + Qi • c, 

a*-^^ = Ki • a^^^ + K 2 • [ad) + b] 0 [a^) + b] + Q 2 • [a^) + b] 0 c. 


(43) 

(44) 


where the dot operator is not merely a tensor contraction, but also involves a time convolution. Eor example, Z being a given 
field: 


(Kl . Z)p(t) dT(Fl)p,q(T - t)Zq(T) , 

a J— 00 


and similarly the higher order contraction rule over the fields Z^ 0 ■ ■ ■ 0 Z"^ is defined as: 
(Kn-Zl0---0Z")p(t)= X] /dTi---/ dTn(k:n)p,qi,...q„(Ti-t,---,Tn-Tn_i)Zl^(Ti)---Z"_^(Tn) . 


(45) 


(46) 


Note that the order of the argument does matter, (i.e. Eq. <461 defines a non-commutative algebra). Note also that the sum of the 
order in each term corresponds to the order of the perturbation. For instance, in Eq. <44< . the second term involves the product 
of two first-order terms, while the first term is a single second-order term. Note finally that the contraction for the Qn involve a 
summation over 5 indices, £, m, a, I', m', (whereas contraction over Kn involves only 3 indices: n, i, m). We illustrate and discuss 
in Fig. 13 through synthetic diagrams the corresponding expansion. (See also Fig. <m< in Appendix O for ari expansion to 
higher order). 
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In appendix we show in Eqs. ID3> - ID4> how to 
rewrite Eq. <44> to order n. 

As for all expansion schemes, the issue of the truncation 
arises. Depending on the physical process investigated, the 
truncation order may vary. Eor instance, it may be legiti¬ 
mate to truncate the perturbation to second-order since the 
second-order is the first-order for which dynamical friction 
is taken into account. 


3.2 Non-linear two-point correlation functions 

Let us now re-address the computation of the two-point cor¬ 
relation function (c/. Section ITU of the response of the halo 
to tidal excitation and infall while accounting for the non- 
linearities described in Section lS.l.ll Eirst, let us reshuffle the 
hierarchy in a format which is best suited for the statistical 
average of the non-linear response. 


3.2.1 Reordering in b and c 

Let us define T(cu, t) = exp{tcvt) the Fourier operator, so that 
T ■ Z and M ■ Z are respechvely the half-Eourier and inverse 
half-Fourier transform of their argument Z. Calling 

■ T, (47) 

Eq. <44> (and its generalization Eq. tD3l l reads like a recur- 
sion: 

a(")=Ri-/C[a(''^^)---,aW,b,c], for n>2, (48) 

where K. sfands formally for some combination of Kn and 
Qn. Nofe fhat Kj accounts for the self-gravity of the halo. If 
the halo is very hot, this self-gravity may be neglected al¬ 
together and Rj 1. If not, we may define Kl = Rj • K;, 
Q' = Rj • Q,', and rewrite the recursive relations Eq. <481 with 
Kj = 0. Eor instance: 

a(i) = Ri • (Ki • b + Qi • c) = K'l ■ b + Qi • c , (49) 

which we can rearrange as: 

a(i) =A^-b + Ac-c (50) 

where Ay = Kj and Ac = Qj. Let us also introduce K" = 
Kj -I-1. Similarly, the contribution of b's and c's fo the sec¬ 
ond order term for a can be expressed as: 

a(2) = Ayy • b (g) b + Acc • c (g c 4- • c g b -f Aye • b g c. (51) 

where 

Ayy = K^oK'/, A„ = K^oQ' -4Q^o[Q',I], (52) 

Acy = k^o[q;,k;'], Afc, = fc^o[K;',Q;] + Q^o[K",i]. 

Here the bracket, [, ] accoimts for the differential composi¬ 
tion, so that, 

• b g b = • (Qi • b) g (K'/ • b) 

= Ri.K2-(Ri-Qi-b)g([l + Ri]-Ki-b). 

In appendix IDl. 21 we also show how to write an equation 
similar to Eq. <511 for the third order contribution and more 
generally for an arbitrary order (see Ea lD8< . 


3.2.2 Non-linear correlators 

We may now complete the calculation of, say, the two-point 
correlation function of the density, C 2 : 

=Y^J^£^{p^^Hxi)P^''^’’\x2)), (53) 

n p=l 

where x,- = (r;, t;), i=l,2. Following Eq. <35< . let us also 
expand the response in density, p, over the basis function 
{pl‘l](r)j.q^ so that 

n qi,q2 

Now, given Eq. l49l and isTl . we may rearrange this equation 
as: 

= £2 E p[‘Ll(ri)p[‘l^l(r2) [cfU£CP + ...] . (54) 

qi,q2 

where is a simple reshuffling of Eq. <29< . i.e: 

= AyxAy • (b g b) 4- AcxAc • (c g c) 4- AyxAc • (b g c) 4- 
AcxAy ■ (c g b) , (55) 

and: 

= (A^i,xA^4-AjxA^i,) • (bgbgb) 4- 

{AccX-Ac 4“ AeXA CC ) • (egege) 4- 
(AycXAc 4- AyxAyc) • (b g c g c) 4- 
[AyyxAc 4- AyxAyc) • (b g b g c) 4- 
(AcyxAy 4- AcxAyy) • (c g b g b) . (56) 

The X operator is non-commutative and guaranties that the 
order is preserved in the dot contraction. Recall that Ay = Kj 
and Ac = Qj, while Ayy, Acc, Acy and Aye are given by 
Eq. <531 (or in terms of the underlying distribution function, 
Fo(I), and the basis function, i/;["l(r) via Eqs IrTi . <421 and 
<4:7< through the definitions of Kj, Qj, K 2 and Q 2 ). If follows 
from Eq. <56< that the non-linear two-point correlation will 
involve at least the three-point correlation of the incoming 
flux and of the external potential. We will see in Section 
that this is a generic consequence of mode coupling. Now the 
three point correlation of fhe incoming flux, c, and the tidal 
field, b may be reexpressed in terms of fhe mean and the two- 
point correlations of thos e fields while relying on W ick's fhe- 
orem, since we showed in lAuhert fe Picho^ ^2^05a^ that these 
fields were approximately Gaussian. Appendix |D2| presents 
formally the generalization of Eqs. <55<-f5^ for the N-point 
correlation function to arbitrary order. 

Equations such as Eq. <441 or its reordered version 
Eq. <51< might look deceivingly simple. One should never¬ 
theless keep in mind that the perturbation theory involves 
an exponentially growing number of terms. This is proba¬ 
bly best realized by looking at diagrams such as Fig. <n2< 
(presented in the appendix) while keeping in mind that each 
straight line represents a triple sum over k = {ki,k 2 ,k^) 
and a time integral (see also Appendix IbI. The prospect of 
achiev ing resummation ( in the spirit of what was achieved 
by e.g. iBernardearJ Il992ll for the gravitational instability of 
the large-scale structures) given the relative complexity of 
the double source expansion is slim. Yet it might be possi¬ 
ble to construct scaling rules (see lErvHl984h l since gravity is 
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also here the driving force. Let us stress once again that the 
perturbative expansion accounts explicitly, within its conver¬ 
gence radius, for all aspects of the non-linear physics taking 
place within the J ?200 sphere. 


3.3 Implication for dynamical friction and tidal stripping 

One of fhe possible assets of this perturbative formulation 
is that the incoming flux may describe a virialized object 
which has a finite extent, and as such will undergo inter¬ 
nal phase mixing reflecting the fact that different points in 
the object will describe different orbits, at different frequen¬ 
cies (see Fig. In the perturbative regime, dynamical fric¬ 
tion will also account for both the overall drag of the object, 
but also its tidal stripping {i.e. the fact that the less bound 
component of fhe object will imdergo a differential more ef¬ 
ficient friction). Specifically, the deflection of perturbed tra¬ 
jectories will correctly describe the balance (or lack thereof) 
between fhe self-gravify of the entering flow and its tendency 
to be torn by the differential gravitational field of the halo 
(which imposes the imperturbed different orbital trajecto¬ 
ries). As such, the flow paradigm implemented in th is paper 
and in I Aubert et alJ i2004h . lAubert & PichonI ^200.5a^ should 
allow for the appropriate level of flexibility in defining whaf 
a structure is and how time-dependent the concept is, within 
the self-gravitating halo (see also Section l3.3.1l . 

Let us briefly discuss how to identify substructures 
within the halo. 



Figure 3. displays qualitatively a bundle of orbits (in their orbital 
plane) which undergo dynamical friction and phase mixing within 
the R 200 radius. As expected, dark matter describing orbits which 
initially are at the same position, but with "slightly" different ini¬ 
tial impact parameters will end up in quite different regions at later 
time. On the right the curve represents a possible angular distribu¬ 
tion of a given entering object (for which the kinematic and angu¬ 
lar spread has been greatly exaggerated). The caustics correspond- 
ing to the successive reboun d of the orbits is clearly visible here 
iFillmore & Goldreichl Il984ll '). Note that the amplitude of the fric¬ 
tion force was ad hoc, and the self-gravity within the bundle was not 
taken into account. 


3.3.1 Substructure counts and distribution 

The identification of subsfructures within a given halo is 
a very promising but difficult topic (see e.g. ISpringel et alJ 
j200lMCill et'n]i2004lAlIbert et alJ<2004 l . Once the bound¬ 
ary flow has been propagated inwards, we have in principle 
access to the full distribution fimction of fhe perturbation as 
a fimction of time. When the field /(v,r, f) is known inside 
1 ^200/ we may attempt to identify collapsed objects and apply 
some form of count in cell statistics in order to characterize 
their spatial distribution as a function of time. This would al¬ 
low us in particular to put aside objects which have been dis¬ 
rupted by tidal stripping or phase mixing (indeed 10 % of the 
mass of the halo is believed to remain in the form of virial¬ 
ized objects, while 90 % is disrupted by the tidal field). Recall 
that the disruption process is in principle well described by 
the perturbative expansion. 

The criterion for the detection of objects must be carried 
while accounting for both the density contrast and the corre¬ 
sponding velocities (see also lArad et alJ i2no4ti . Indeed, we 
do not wish to identify as objects local overdensities which 
may just correspond to caustics or local wave reinforcement. 
Here we are interested in the temporal coherence of objecfs. 

For this purpose, we may coarse-grain the perturbed 
distribution function both in position and velocity, with some 
given smoothing function, W(r/72s, v/Vs), and then apply 
some thresholding (W o / > where o stands for convo¬ 
lution) on the amplitude of the distribution function, defin¬ 
ing a set of coimex regions. For each of these regions, we 
may then compute the energy of the corresponding clump. 
If it is negative, the clump will be labeled as bound for the 


corresponding threshold (/min)/ and coarse-graining param¬ 
eters (72s, Vs). Note that since the response only involves the 
perturbed density, one need not subtract the mean potential, 
(which is quite a difficult task in general). 

Once the bound regions are identified, we may compute 
the corresponding mass and assign it to the bottom of the lo¬ 
cal potential well. This procedure may be applied for a range 
of threshold values, and standard statistical tools for discrete 
sources but in spherical geometry. We may in particular con¬ 
struct in this manner the mass function of satellites as a func¬ 
tion of radius, or, say the two-point correlation function ver¬ 
sus mass and cosmic time. Both issues are subjects of strong 
discussions when addressed through standard N-body sim¬ 
ulations. 

This time-dependent identification of virialized objects 
is useful because of biasing, i.e. the fact that most observa¬ 
tional tracers will only be sensitive to the more massive tail 
of the mass function of virialized objects. 

Conversely, we may want to label regions which match 
the thresholding but not the requirem ent on binding energy, 
i.e. ide ntify caustics, cusps, and shells iFillmore & Goldreichl 
il984h L We may then characterize statistically the mean dis¬ 
tance between the apoapses (see Fig. j3), which will in gen¬ 
eral depend on 72s /Vs and but also on the imderly- 
ing equilibrium, via F(I) and on the statistical properties of 
c fhrough, say the distribution of impact parameters. Note 
in closing that the competing effects of phase mixing, tidal 
stripping and dynamical friction all assume fhaf the imder- 
lying basis function reaches sufficiently high spatial frequen¬ 
cies to resolve these phenomena. In practice, since the pro¬ 
jection of the response (both linearly and non linearly) is 






















Dynamical flows through Dark Matter Haloes 11 


achieved over a basis which has a truncation frequency, fmax/ 
there is a finite time scale, T^ax imaxi (<^) above which 
phase mixing would induce winding at unresolved scales 
(here (w) represents the typical frequency of the dark mat¬ 
ter in that region). Since the dynamical time is shorter in the 
inner region of the galaxy, such a threshold is going to be 
reached there first. Beyond this critical time, the dynamics is 
inaccurately modeled for the corresponding clump. This is¬ 
sue will be important for the non linear coupling of clumps, 
since substructures entering the halo at later times will be 
dragged by streamers which are beyond the accuracy thresh- 
old^ 


4 QUASI-LINEAR EVOLUTION OE HALO PROEILE 

The previous sections dealt with the halo polarization while 
considering that perturbations were transients. In practice, 
a halo undergoes recursive excitations from its environment 
that will induce departures from its equilibrium state so that 
it won't remain static. In Appendix^we derive a quasi lin¬ 
ear formalism for the collisionless open Boltzmann equa¬ 
tion in order to tak e this effec t in ac c ount. This follow s in 
essence the work ofjj ^inbergl il99,'^i . I Weinberg i2nniah . or 
& Bertschinge j i2004h . though the derivation differs. We 
introduce here an explicit expression, valid at low redshift, 
for the source of stochastic "noise". We account explicitly for 
the correl ation induced by the ent ering material (as charac¬ 
terized bv lAuhert fe Pichorl i2fln,SalB rather than rely on some 
ad hoc assumption on its nature. We also account consistently 
for the mean secular infall which adiabafically restructures 
the mean profile. 


4.1 Context and derivation 

iGilbertI <197011 gives a very elegant derivation from first 
principles of the secular equation based on a 1/N (N 
being the number of particles in the system) expan- 
sion of the collisiona l relax a tion equat i ons p re sented by 
iBogol yubov & Gurov! <1947^. IWeinberd <199.ift . IWeinberp 
<2001ah and iMa & Bertschinged <20n4h relv on the same ex¬ 
pansion scheme to derive their kinetic equation for the mean 
hal o profile. 

IWeinber3 <199.'jl focuses on the secular collective relax¬ 
ation of a system induced by the finite number of particles 
within a multi-periodic uniform medium, hence transpos¬ 
ing to co llisionless stel l ar dynamics the derivation of Lenard- 
Balescu <Lenardl <196lh . l^lescul Jl96.'jl 'l applied originally to 
plasma physics in order to describe the secular convergence 
of s uch systems towar ds thermalisation. 

IWeinber^ <2001ah derives a similar result for the spher¬ 
ical halo in angle and action va riables, while relying on 
the Kramers-Moyal <Riskerj lI^l expansion, which corre¬ 
sponds to a Markovian description based on the transition 
probability of a change in action induced by the interaction 
with a dressed particle cloud. His Fokker-Planck coefficients 
differ slightly from Eqs. <Cil.3> - <rd^ given in the Appendix 


* these limitations are also clearly encountered in classical N-body 
simulations 


in that the spectral properties of (fen^n') are postulated in his 
case , while c = 0. 

iMa & Bertsching^ <2n04h construct a Fokker-Planck 
equation for the mean profile of a halo in a cosmic environ¬ 
ment while relying on the constrained random field of peaks 
in the standard cosmological model to derive the drift and 
diffusion coefficients from first principles. Their derivation 
is dynamically accurate to second-order in the perturbation 
theory (in position-velocity space) and relate the kinetic coef¬ 
ficients to the properties of the underlying linear power spec¬ 
trum. In contrast to the theory presented here, their kinetic 
equation describes the very early phase of halo formation, 
whereas we focus here on the quasi-linear evolution (in an¬ 
gle action space) of fully relaxed equilibria at low redshift. 

In Appendix 10 we account explicitly for the na- 
ture of the perturbation 's power spectrum as defined in 
lAubert & Pichon| <2n0.5ah and present an explicit derivation 
for the Fokker-Planck equation obeyed on secular time scales 
by the distribution function in angle action. It is natural to 
use these variables to describe a relaxed collisionless halo 
since they allow to split the dynamics into a secular (phase 
averaged) and a fluctuating part. 

Even though individual dark matter particles obey a col¬ 
lisionless dynamics, the phase average ("ensemble average") 
distribution for the open system satisfies a collisional kinetic 
equation where the dumpiness of the o pen medium breaks 
the m ean field approximation (see also iMa fe Bertsching^ 
<2n04h '). Indeed, individually, clumps and tidal remnants de¬ 
flect the actions of fhe underlying distribution in a stochastic 
(but correlated) manner, so that in the mean ensemble sense, 
the coarse-grained distribution {i.e. the distribution averaged 
over the angles) obeys a collisional diffusion of the Fokker- 
Planck type. In this formulation, the graininess of the system 
(as defined by fhe second-order closure of fhe BBGKY hier¬ 
archy of the N-point distribution) corresponds to the mean 
number of clumps expected in the halo, while the detailed 
(kinetic and angular) power spectrum of the gravitational 
fluctuations is given by the cosmogony. 

It is usual in plasma physics to take a two-time scale 
approach to the Boltzmann equation. The short time scale 
describes the system's dynamics on the dynamical (orbital) 
time scale, while the longer time scale corresponds to the sec¬ 
ular evolution. The Action-Angle variables are best suited 
here. This time scale separation procedure leads to the fol¬ 
lowing system of equations: 


dt 9w 3w 91 


9T \ 9w 9w 



dlpe 9F 

9w 91 


dtp dtp'^ ’ 9/ 

91 91 9w 


(57) 
+S„ (58) 


In Eq. <571 and <58t Se and Se stand for the perturbative and 
secular adverted source terms, while / stands for the fluc¬ 
tuating distribution and F stands for the secular distribution 
function (see Appendix O for defails). The bracket around 
the quadratic terms stands for a time average over a sec¬ 
ular time, T which is long compared to dynamical time, t 
(taken by a dark matter particle to describe its orbit). If we 
fix F(I, T), Eq. <571 corresponds exactly to Eq. whose so¬ 
lution was described in Section lTTl This formal solution may 
then be injected in the quadratic terms of Eq. <581 . Following 
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Dark matter deflection on polarised 
correlated clumps in position velocity (left) 
and angle-action (right) space 


Figure 4. Left panel schematic representation of the successive deflection of a given orbit on correlated clumps within the halo in position space. 
Each clump is represented with its polarisation cloud. The gray scale coding in the cloud reflects their spatial and temporal correlation within 
the clumps. During the deflections, the orbital parameters change (though the individual change is here grossly exaggerated). Right panel: 
the same orbit as viewed in angle-action variables. The dynamics in these variable is straightforward (it corresponds to straight lines obeying 
w = c<;(I)f + Wfl) and the diffusion process resembles Brownian motion obeying a Langevin equation the particle receiving a random kick at 
each deflection, represented by a change in colour which reflects the fact that the "collision" is instantaneous in constrast to the time interval 
separating two collisions 


this route, we show in AppendixOhow to rearrange Eq. i58l 
as a Fokker Planck equation: 

— = (Do(I)) - (Di(I)) . — - (D 2 (I)) : (59) 

where Dg, Dj and D 2 are given by Eqs. IC13l - td^ . while 
: stand for the total contraction. Note that Eq. <591 . in con¬ 
trast to Eq. <58< refers this time to the driving equation for 
the mean halo profile since we invoqued ergodicity to replace 
time averages by ensemble averages (see AppendixOfor de¬ 
tails). The Dg term enters here because the halo is an open 
system, which may receive or lose mass. The drift term with 
the factor Di accounts for the dynamical friction induced by 
the polarization cloud aroimd the tidal remnants; the diffu¬ 
sion term with the factor D 2 arises because of the fluctuations 
in the potential (both tidal and associated with the infalling 
dark matter) induced by the clumps. The diffusion term will 
in general induce a spreading of the energy distribution by 
accelerating some orbits to higher energies while decelerat¬ 
ing some other ones. The polarization cloud will in general 
induce a drag on the clumps, represented by the Dj term. 
Note that the former should be independent of the mass of 
the clump (since the energy is exchanged via the mean field) 
while the latter will not (since more massive clump polarize 
the medium more). From the point of view of the entering 
dark matter, the net effect is therefore a segregation process 
in which the more massive clumps fall in (as discussed by 

ICilbertlJl97Cl) ). _ _ 

According t o lRiskerJ il989ll . the corresponding Langevin 
iLangevinI h908ll ) equation reads (when the source term, Dg 
is omitted) 

| = Ai(I) + A2 (I)-?(t). (60) 

Here Ai(I) and A 2 (I) are given in terms of Di(I) and D 2 (I) 
by 


A 2 = , and Ai = Di - , 

where [D 2 ]^^^ stands for the square root of the matrix D 2 
which is computed via diagonalization, provided the eigen¬ 
values are positive. The 3D random field, f(f), should have 
spectral properties which reflect the stochastic properties of 
b and c. The probability distribution of the solution to the 
stochastic equation, Eq. <60< , obeys the Fokker Planck Equa¬ 
tion, Eq. In this form, the effect of diffusion on the de¬ 
parture from phase mixed equilibrium is easily interpreted. 

4.2 Prospects for universal halo profiles 

As has been suggested and il lustrated bv lWeinberd ^2001a^ 
and iMa & BertschingeJ i2004h . it would be very worthwhile 
to use Eq. and predict the asymptotic dark matter profile, 
(and, say the cosmic evolution of the concentration parame¬ 
ter) which will be shaped in part by encounters and interlop¬ 
ers. 

Note that the diffusion coefficients, D; are relatively 
straightforward to compute for a given halo model, f (I) but 
Eq. <591 corresponds to an evolution equation for F(I) and 
will in practice require re evaluating the coefficients for dif¬ 
ferent values of F. 

Let us now draw constraints on the stationary solutions 
of Eq. <591 . Again (following Section|3, this may be done in 
one of two ways: take Dg, Dj and D 2 as given function of 
the actions, and deduce what equation F should obey from 
requiring that Eq. <59< has a stationary s olution, (this is the 
route first explored hv lWeinhergI i2001hh ): or, if we assume 
that a given model, say a universal profile, should corre¬ 
spond to the asymptotic solution of Eq. <591 . we may find the 
relationship relating the corresponding asymptotic D,- coeffi¬ 
cients. 

For simplicity, let us illustrate this second point while 
neglecting here the fact that the diffusion coefficients depend 
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on the distribution function, and restricting ourselves briefly 
to an isotropic distribution, F{E, T). Calling: 

: / ((D 2 ) : 07® w), (61) 


H{E)= (Di).a7 + (D 2 ) 


and Q(£) 


(Dp) 

(D 2 ) : 07 ® 07 ' 


the stationary solution {dF/dT 
mally: 


(62) 

0) to Eq. reads for- 


F{E)= f exp — f H[ei]dei 

Jo . J ^0 


QM exp 


Jei 




(63) 


This distribution function should satisfy the self-consistency 
requirement that 

fO 

p{r)=2V2 F{E)^/ETlpdE, V^xp{r) = inGpir). (64) 
J — tp 

Imposing that F(£) obeys Eqs. i^- <64> yields a non-linear 
integral equation for the D;, i.e. a (admittedly indirect) con¬ 
straint on the angular correlation of the external field. Sec- 
tionj^describes other means of c onstraining the po wer spec¬ 
trum of the infalling dark matter. lWeinhergH2n01 bh foimd it¬ 
eratively the corresponding solution while making some as¬ 
sumptions on the spectral properties of b in the regime where 
c = 0 . In the light of his investigation, he concluded that the 
tidal excitation drives the halo towards a less steep profile. 
It will be interesting to explore this venue with a realistic ac¬ 
counting of the source of infall. The sethng here would be 
that the satellite problem and the cusp problem of dark mat¬ 
ter haloes might be the two sides of the same coin, so that the 
evolution towards a universal profile might be triggered by 
the actual infall of substructures. 

Let us now return to the perturbative dynamics de¬ 
scribed in Section 12113 and explore its implications for galax¬ 
ies. 


5 APPLICATIONS: HALO POLARIZATION, DISK 
DYNAMICS AND INVERSION 

lAubert & Pichon| i200.'5allJ) provided a detailed statistical de¬ 
scription of how dark matter falls onto a L* galactic halo: 
how much mass is accreted as a function of time, how is it 
accreted i.e. in what form, with what velocity distribution, 
along which direction, and for how long? Putting the theory 
described here and the tabulated measurements from that 
paper together, allows us to address globally, and coherently 
dynamical issues on galactic scales in a statistically represen¬ 
tative mariner. With the help of the theory presented in Sec- 
tion|2||3 we are now in a position to ask ourselves: what are 
the expected features of a halo/galaxy induced by their cos¬ 
mic environment. Specifically, we may now "simply" prop¬ 
agate the cosmological framework and its statistics to ob¬ 
servables (describing the departure from spherical symme¬ 
try/stationarity) on galachc scales. On these scales, the realm 
of astrophysical applications for the perturbative open solu¬ 
tion of the Poisson-Boltzmann equations is extremely wide. 
It is clearly beyond the scope of this paper to attempt an 


exhaustive inventory. Rather, we shall here focus on a few 
specific issues, for which we show how the open perturba¬ 
tive framework improves our understanding, and allows for 
a statistical investigation. 

In particular, we shall restrict ourselves to settings where 
the detailed geometry of the infall matters, since the theory 
described above does account for the configuration and the 
time lag involved in the accretion on top of L* galaxies. 

Recall that the purpose of the statistical propagation is 
threefold: (i) constrain the properties of the infall on the basis 
of the observed distribution for the properties of galaxies and 
their environment; (ii) predict some of the statistical proper¬ 
ties of galaxies which are not directly observable, while re¬ 
lying on the properties of the infall, (iii) weigh the relative 
importance of the intrinsic properties of the disk+halo com¬ 
pared to the strength of the environment. 

We will distinguish three classes of problems; first we 
will describe how to transpose to galactic scales fSection BTl 
the classical probes used in cosmology to trace the large-scale 
structures. We will then explore in Section lSi^ the implication 
for the properties of external galaxies, and in Section I531 for 
the structures within the Milky Way halo. Finally, we will 
elaborate in Section ISTl on the prospect of inverting the up¬ 
coming data sets for the past history of our own Galaxy and 
for field galaxies in the local group. 


5.1 Cosmic probes in the neighbourhood of galaxies: 

1^200 /10 < R < R 200 

A series of observational probes of the statistical properties 
of the density field have been devised over the years, such 
as weak lensing, galaxy counts, the SZ effect. X-ray or 7 - 
ray emissivity maps. In the light of large galactic surveys 
which are available today, it becomes quite desirable to apply 
these probes in the neighborhood of galactic haloes in order 
to study the dark matter distribution within the R 200 radius. 
Some of these tracers are only sensihve to the baryon den¬ 
sity, which need not trace directly the dark matter density. 
In this section, we will systematically assume for simplicity 
a simple biasing, though this assumption may be lifted (at 
the expense of extra non-linearihes, see Section lE51 provided 
the biasing law is known (i.e. the observables are assumed to 
scale like the dark matter density, or some power of it); we 
refer to Section 13.3.11 for a brief discussion of thresholding, 
which is boimd to be important in practice. 

The calculation described in the previous sections, 
together with the statistical measurements described in 
lAu^bert et alJ i2004h : lAubert & PichonI i2005allJl should allow 
us to make statistical predictions about observables which 
may be expressed in terms of the distributions of clumps 
within the galactic haloes, either via their gravitational po¬ 
tential, their projected density or even their velocity distri¬ 
butions {e.g. Galactic streams). 

We will consider in turns observable which may be ap¬ 
proximated as linear fimctions of the perturbed fields, either 
in projected coordinates on the plane of the sky, or as seen by 
an observer at the galactic center. We will also consider ob¬ 
servables which involve quadratic functions of this field {e.g. 
the square of the electron density), or even more non-linear 
functions of the dark matter distribution within the virial ra- 
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dius (such as the locus of virialized clumps, which dissolve 
at a function of time). We will in particular build the two- 
point statistics for these observables, since the mean of the 
perturber is often zero by construction. Finally we will also 
consider metals lines in absorptions systems, which involve 
the cross-correlation of fhe density and the velocity fields. 
Note that all these measurements could in principle be car¬ 
ried as a function of redshift, or a function of fhe mass of 
the halo, or while varying the anisotropy of the equilibrium 
for the halo (by varying F(I) in e.g. K in Eq. im v Note fi¬ 
nally that some tracers correspond to the scales of clusters, 
and we will assume here that the measurements presented in 
lAubert & Pichon| i200.5ah could be reproduced for these ob¬ 
jects (whereas the theory described here is scale independent 
provided the system is dynamically relaxed and spherical). 

In this section, we will focus on a couple of probes which 
are supposed to scale linearly with the dark matter density 
in the main text (weak lensing, SZ effect), and postpone to 
appendixl^a presentation of other probes (X-ray emissivity, 
dark matter disintegration, metal lines in absorption spec¬ 
tra). 

Note that all probes described below are a departure 
from the mean profile of galactic haloes (just as cosmic per¬ 
turbation theory describes the growth of structure as a depar¬ 
ture form the mean density/expansion of the universe) and 
as such, assume that we have a good understanding of this 
profile. This will undoubtedly turn out to be a serious ob¬ 
servational constraint when attempting to ensemble average 
galaxies of various size and properties. 


5.1.1 Weak lensing in stacked haloes 

Weak lensing corresponds to the deflection of light emitted 
from background galaxies by the gravitational potential of 
structures between those galaxies and the observer. It has re¬ 
cently been used quite successfully to constrain the statistical 
distribution of the large-scale structures. On smaller scales, 
the effects of substructures in haloes o n the lensing mea- 
surments have been demonstrat ed by e.g. lDalal & KochanekI 
<2002h . lKochanek & Dala]l<2004 . 

In the weak lensing regime iPeacockI iigggh l. the re¬ 
lationship between the observed convergence and the un¬ 
derlying projected dark matter profile is approximated to 
be linear. Hence we may straightforwardly propagate our 
statistical predictions for the clumpy dark matter distribu¬ 
tion around a dark matter halo (or within the neighbor¬ 
hood of clusters of galaxies provided s ome readjustment o f 
the theoretical predictio ns described in lAubert et alJ i20n4l : 
lAubert & Pichon| i2005ah on these larger scales). 

The cumulative deflection angle, x(9,w) = Sx/rpzv) by 
which light is deflected is given by 


the mean ensemble average convergence of the rescaled halo, 
(k( 0 )), reads: 

(K(e)) = ^ j (67) 

Now recall that (c/. Eq. where tpl"! is given by Eqs. <B2I - 
tB3l in Appendix lB2l ') 

Jipjr) = ^ hence SK{d) ='^ar,P''\d), ( 68 ) 

n n 


where 


KM(e) = \ f (69) 

J ri,{w) 

It follows fhat the correlation function of the relative conver¬ 
gence obeys 


(<5k(0)<5k(0')) 

( k ( 0))2 


1 

W¥ 


^{ar,a^,)P'^\e)P'^\0'). 

n,n' 


(70) 


Hence, the statistical properties of fhe relative convergence 
will depend on the statistical distribution of fhe clumps of the 
halo through the {an} coefficients which are given in Eq. 
in terms of bn and Cn. 

In practice one has to devise an observational strategy, 
given the expected size of fhe caustics of subclumps within 
haloes of galaxies or clusfers, the number of background 
sources, and the expected number of foreground objects (i.e. 
galaxies or clusters). 

Einally, it is believed that one in a hundred large ellip¬ 
ticals on the sky should undergo strong lensing. In the long 
run, the statistical properties of such a non-linear signal will 
be worth investigating within the framework described in 
this paper (following the non-linear steps described in say, 
SectionlEsl. 


5.1.2 Thermal S-Z effect of stacked haloes 

When the photons of fhe cosmic microwave backgroimd en¬ 
ter the hot dense gaz within the clusters and galachc haloes, 
they interact with the electrons of fhe gaz. The diffusion 
process transfers fhe energy of the photons to the electrons 
which in turn reemit this energy at a higher frequency. The 
corresponding spectral redistribution induces a local temper¬ 
ature decrement seen in the temperature map of the clus- 
ters . known as th e thermal Sunyaev-Zeldovich effect (see 
e.g. lPeacockI il999h l. The temperature decrement (at low fre¬ 
quency) reads as a function of the distance to the cluster cen¬ 
ter, R : 


a(0,ry) =— ^^57^xp{rpzv')d,w'), (65) 

c2 J r^{w) 

where rj- is the angular comoving distance, 
dw = dr/ Vl — kr^ {k = 0, ±1) and x the transverse co¬ 
moving distance. Defining the convergence, k( 0) by 

K(0) = lvfl-a(0), (66) 


AT(R) 

Tcmb 


' mgC^ 


J dzwe(z, RjTgjz, R), 


(71) 


where m^, We and Tg are respectively the mass, the numerical 
density and the temperature of the electrons, while Oj is the 
Thomson scattering section (6.65 x 10^^®cm^), c the speed of 
light, fcj Boltzmann's constant, and T^mb is the Cosmic mi¬ 
crowave backgroimd temperature. 
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Let us assume that the variation in temperature is small 
compared to the variation of the electron number density^. 
Let us also assume that the electron density is proportional 
to the dark matter density (constant biasing) as mentioned 
above. Let us define the departure from the cosmic average 
for the profile as: 

<5AT(R) = AT(R) - (Ar)(R). (72) 


The relative fluctuahon of the temperature decrement reads: 


1 


-(MT(R)MT(R')) =: 


^{anan') 


(AT)2(R)^ V ^ V 

J dz J dz'p[”l(R,z)p["^(R',z^), (73) 


where Snfw is the mean rescaled projected dark matter 
mass profile. Note that the double integral in Eq. <731 is car¬ 
ried over known functions and is just a geometric factor 
which will depend on R, AR, and A© only. Again, the knowl¬ 
edge of the statistics of the {fln} (which in turn only depend 
on the equilibrium, Fg, and the statistics of fen and Cn at R 20 O/ 
see Eq. or 03) therefore allows us to predict the sta¬ 
tistical properties of the relative fluctuations in the tempera¬ 
ture decrement. ALMA will soon provide detailed SZ maps 
of clusters for which it should be possible to apply these tech¬ 
niques. 

Let us note in passing that the kinetic S-Z effect of 
stacked haloes may also be investigated following the same 
route 


=-2^ [ dzne{z,RM{z,R). (74) 

LcmB nieC^ J 

In closing, let us note that maps of SZ effects within 
our own Galaxy will be available with the upcoming Planck 
satellite, and will provide statistical information on the 
small-scale distribution of local clumps. Recall finally that 
appendix presents other statistical probes of the outer 
structures found in galacfic haloes (X-ray emissivity, dark 
matter disintegration, metal lines in absorption spectra). 
Most of these probes could be used to say, probe the shape 
of the density profile in the outer parts of galaxies, or the bi¬ 
asing law relating dark matter to stars or gas. 


5.2 Galactic structure: R < R 200 /IO 

In the previous section, we investigated the dynamical con¬ 
sequences of the cosmic infall in the outer region of the halo. 
Let us now turn to the regions of the halo where we expect 
to find the galaxies themselves. At lower redshift, the galax¬ 
ies essentially come in two flavours, ellipticals and spirals. 
The response of ellipticals should follow closely that of the 
dark matter halo since both components are hot enough not 
to undergo gravitational instabilities. In effect, describing an 
embedded "spherical" elliptical galaxy within a dark mat¬ 
ter halo amounts to changing the distribution fimction to ac¬ 
count for the presence of the elliptical and its possibly dis¬ 
tinct kinematics. Note that, as mentioned before the above 

^ Note that we may lift this assumption at the cost of nonlineari- 
ties, provided we may rely on an equation of state to relate it to the 
underlying density 


theory could be amended to acco unt perturbatively for th e 
possible triaxiality of the elliptical iBirmev & Spergeilh984h L 

Eor a disk or a very flattened spheroid, the situation be¬ 
comes quite different. The cooler disk is likely to be either 
drawn beyond its stability threshold by the perturbation, or 
will respond much more strongly to the perturber than its 
dark halo. Hence we need to model the disk component dif¬ 
ferently. The proto-galactic environment is likely to be ex¬ 
tremely noisy, particularly in outer regions, so that the halo 
may perturb the disk by transmitting numerous disturbances 
into the inner galaxy. Moreover, the inner halo may continue 
to oscillate as it settles after the coalescence of advected ob¬ 
jects. Halo oscillations may easily perturb the disk through 
the time-dependent gravitational potential. Conversely, the 
structural integrity of observed disks set limits on the degree 
of disequilibrium in the proto-galactic halo. 

Within the realm of features found in galactic disks, a 
fraction are known to be the result of instabilities {e.g. galactic 
bars), while others have been shown to correspond to tran¬ 
sients (e.g. galactic warps). 

With the advent of modern systematic surveys, it is pos¬ 
sible to construct distributions corresponding to, say, the 
fraction of spirals which fall within some Hubble type, or 
the fraction of warps whose inclination is larger than some 
angle. On the disk scale, we may construct the PDE of, say, 
the pitch angle of dynamically induced spirals, or the PDF of 
the extent of the bar, its amplitude, or less directly observed 
the PDF of pattern speeds. Some of these processes depend 
crucially on gas physics and will not be addressed here. 


5.2.1 Pitch angle distribution for spirals 


For stellar disks, the stars obey formally the same equation 
as Eqs. l^-fiTl. but this time the modes may be unstable, 
and sometimes the disk cannot be treated in isolation from 
the live halo in which it is embedded. On the other hand, 
it is often well approximated as an infinitely thin structure; 
such a 2D system becomes integrable again with two actions, 
dj = d/rdLz p. Here f is the radial action of the stars in the 
plane of the dis k, and Lzr, is the mo mentum of the stars in the 
disk. Following I Weinberd ^1998a^ and adding some source 
of infall at R 20 O/ we may describe the coupled system disk + 
halo in the complex plane as 


/ ao A / Kdo Vr “D V 0 
KM J y Koh Khh j y SH+b ; Mv 0 

where Khh is given by Eq. <D15> . while 


(KDD)p,q 



k • Wp — a; 



0 

c 


and a similar expre ssion involving for the 

cross term, Khh- See lPichon & CannonI il997h for details rel¬ 
ative to the disk. Here ao and an are the coefficient of the 
expansion for the disk and the halo respectively, F„ is the dis¬ 
tribution function of stellar stars within the disk, {t/jD'!"!] (r)}q 
the potential basis function over which the disk response is 
projected, and the angular frequencies of the stars in the 
disk. 

Let us first assume that the unperturbed disk is stable. 
Solving the coupled equation for [Sh, Sh] yields (after half 
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inverse Fourier transform) the temporal evolution of the spi¬ 
ral response as a function of time for a given tidal field, b(t) 
and a given infall history, c(t). The pitch angle, I of the spiral, 

definedby tan(J) = I/tt [ d0dlog72./d6 (where7?.(6) cor- 

Jo 

responds to the crest of the spiral wave), is a non-linear func¬ 
tion of ao, which we may write formally as ir[ao]. Hence we 
may ask ourselves what it's cosmic mean, (X[aD]) is, given 
that aD(m) obeys 

= Khd • [b -h Q • c] 

Note that X[ao] will depend on the statistical properties 
of b, c and also on the distribution function for the halo, f (I), 
and the distribution function for the disk, (J). More gener¬ 
ally, we may in this manner construct the full PDF of the pitch 
angle, as a function of say, cosmic time, (or relative mass in 
the disk or ... ) following the same route as sketched in Sec- 
tion lF3l 

If the disk is intrinsically unstable, we must then add 
to the driven response described above the unstable modes. 
The amplitude of the response will then depend on exactly 
when each unstable mode has been exited. Such a prescrip¬ 
tion is beyond the scope of this paper, but could be addressed 
statistically through the description of a phase transition. 


Det 


1 y; Kdd 

Kdh 


K; 


DH 


K, 


HH 


5.2.2 Warp excitation 

As mentioned ear lier ihopez-Corredoira et alJ i20n2ll . 
Ihang & BinnevI <19991) 1. warps are intrinsically stable modes 
of thin disks which respond to their environment. The action 
of the torque applied on the disk of a galaxy is different for 
different angular and radial positions of the perturbation. 
The warp's orientation and its amplitude are functions of 
the external potential. 

The work done by the presence of perturbations on the 
stellar system is 

^ = - y drV {ip + xp'^) ■ pv = - J dt{xp + f)\/{pv), (75) 

where xp + xp‘^ is the total potential perturbation (self-response 
+ external component). 

Using Eq. t20l and Eq. j75j 

(^) = -X] /dr([«n'(t) + VW]Vi/-„(r).x 

n,n' 

f dT]<fp],„(r, T - t) [fln(T) + bn(T)] + Q[2],„(r,T- t)Cn(T)). 

— CO 

The power spectra of potential fluctuations drive the en¬ 
ergy rate of change through the cross-correlation between the 
source and the potential. 

Note in closing that the framework described in this pa¬ 
per should allow us in the future to address the possibility of 
warps induced by the accretion of gas. 


5.3 Substructures in our own Galactic halo 

Let us now turn to the Milky Way. Our knowledge of the 
structure of its halo has increased dramatically in the course 
of the last decade with the advent of systematic imaging and 


spectroscopic surveys (e.g. SDSS, 2dE), both in the optical and 
at longer wavelengths {e.g. 2MASS), and this observational 
investigation will undoubtedly continue with efforts such as 
RAVE, or the upcoming launch of GAIA. This has led to the 
discovery of quite a few substructures within our halo, both 
in projection on the plane of the sky (tidal tails) as star counts 
but also via kinematical features (streams). The extent of the 
upcoming systematic stellar surveys will allow for a system¬ 
atic analysis of the dynamical properties of Galactic substruc¬ 
tures. 


5.3.1 Extent of tidal tails and streams in proper motion & 
galactic coordinates 


The number of stars, dN, in the solid angle defined 
by the Galactic longitudes and latitudes {£, b) =£ (within 
dM(sin&)), with pro per motions (uf, Uh ) —ft (within df(j,df(^) 
at time t is given by jpirhnn et alJ i2nn2h 'l: 

dN = Ax{fi,£,l) = ^JJ<iurPdif{i,u,t)'^ dfidl, (76) 

The variables r, u are the vector position and velocity coordi¬ 
nates {Uy, III, Mb) in phase space relative to the Local Standard 
of Rest, while r = (R, <J>, z) and v = (e>r, v^, Vz) are those rel¬ 
ative to the Galactic centre. In particular, the radius r (within 
dr) corresponds to the distance along the line of sight in the 
direction given by the Galactic longitudes and latitudes (£, b) 
(within the solid angle d^ cos(i’)db). 

These velocities are given as a function of the velocities 
measured in the frame of the sun by 

1 




R 


(rQ sin(5) sin(t')Mf, — Tq cos(b) sin{£)ur— 


rQ cos(£)ui + rcos(b) (ui — sin(£)MQ) -t- (rQ -f rcos(b) cos(f)) Vq) , 
wr = ^ {(r cos(b) — rQ cos(£)) sin(b) Uj, — rQ sin(£) ui— 

cos(b) (r cos(b) — rQ cos{£)) iir+ 
rQ Uq — r cos(fe) cos(f') Mq -b r cos(b) sin(.f) Uq} , 
sm{b) Ur + cos{b) Uf, + Wq , (77) 


Vz = 

where 
= tan“ 


1 /r cos(fe) sin(£) Cq — i cos(b) cos(£) 


R 


R 


R = Y/?^^^-2rQ77os(byTos(£y+^^Tos(^^, and 
z = rsin(fc). (78) 


R measures the projected distance (in the meridional plane) 
to the Galachc centre, O the angle in the meridional plane be¬ 
tween the star and the Galachc centre, while z is the height 
of the star. Here Uq, Vq, Wq and Tq are respectively the com¬ 
ponents of the Sun's velocity and its distance to the Galactic 
centre. 

Recall that Eq. |7| together with Eq. <19> provides us with 
the full phase space distribution of the infall as a fimction of 
the actions of the unperturbed halo. Let us call /n, the phase 
space basis defined by 

/n(r, V, t) = X^expjik • cut -I- ck • w) v zk • (79) 

so that the perturbation at time t and position (r, v) reads 
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dT/n(r,V,T- f)fln(T) . 


(80) 


We may now seek the characteristic signature in observed 
phase space (today i.e. at t = 0 ), of a given perturbation. 


= ^JdTanir)JdMjdiirfn{tl£,r],'v[Tfi,iir],T), (81) 


where v[r^, Mr] is given by Eq. <771 . and r[^,r] is given by 
Eq. <781 . In particular, we may compute the autocorrelation 
of the kinematic count defined by 

= {A{£ + M,fi + Afi)A{£,fi)). (82) 


It involves an integral over time of the autocorrelation of the 
coefficients, (fln(T)fln'(TO) 


^dTdT'(fln(T)fln'(T')) JJ T^drr'^dr' JJ dtiydUr 

fn{t[£, r],v[ifi, Mr], t) fn'{t[£ + M, r], v(r'^ + I'Afi, u'], t'). 


Recall that (fln(T)fln'(T')) can be reexpressed in terms of 
the coefficients of (b • b*^) and ( c • c*^) and (b • c*^) via 
Eq. Og. The width of the correlation, C^(A^, Afi), both in ve¬ 
locity space and in position space accounts for the expected 
cosmic size of structures within the Galactic halo. 


5.3.2 Angular extend of tidal tails 

The marginal distribution over proper motions of Eq. <761 
yields the projection on the sky of the perturbation: 


A{l,t) = JJ A{£,fi,t) dfi =JJJ duMdTf{T,u,t)dfi, (83) 

which can be derived from Eq. <811 but is also found directly 
via integration over the density as 

A{£,t) = YlMt) [ |5n(r,^)r^dr, given 

n 

Pn (r, £)^pn{R{i,£), ‘bjr, £), z(r,£)), (84) 


where pn is given by Eq. and e.g. R(r, £) is given by 
Eq. <781 . Note the generic difference between Eq. <81 < and 
<841 : the former involves the explicit cumulative knowledge 
of fln(T) for all T since it involves a kinematical (inertial) 
quantity, fi, while the latter only require the knowledge of 
the current fln(f)- This difference is weaker than it seems in 
practice, since self-gravity implies that fln(f) depends in turn 
on the previous fln(T) via Eq. <18< . The corresponding angu¬ 
lar correlation reads 


CA(Af) = {A{£ + A£)A{£)) = 

X](fln(f)fln'(0) ff didi' r'Vp„[i,£]pn[i',£ +Al]. (85) 
n,n' 

The FWHM of the correlation defined by Eq. <851 corre¬ 
sponds to the"cosmic" width of tidal stream projected on the 
sky. 


properties of galaxies back in time and constrain the past in¬ 
fall and the tidal field on a given dark matter halo. In short, 
the idea is to notice that the perturbation theory provides an 
explicit relationship between the response and the excitation 
of the inner halo which we can tackle as an integral equa¬ 
tion for the source. ® Let us present first the inversion for our 
Milky Way (Section [5.4.1< . and discuss briefly extra galactic 
stellar streams. 


5.4.1 The Galactic inverse problem 

Let us rewrite formally Eq. <81 < as A(£,p) = AH ■ a, where 
the dot product accounts for both the summation over n and 
the integration over t (c/. Eq. <451 1. Let us assume that we 
have access to kinematic star counts, i.e. to a set of measure¬ 
ments {A,- = A{£i,fii)}i^„. We want to minimize 

^" = E(^<-4''''J^i'(K-b + Q-c))", (86) 

i 

subject to some penalty function. Recall that Rj is given by 
Eq. <471 and accounts for the self-gravity of the halo. Let us 
formally rewrite again • Ri • (Kj • b + Qi • c) = M ■ b, with 

b = ]b, c]. Let us also write A = (A;);^„. The solution to the 
linear minimization, Eq. <86< reads 

b = M[^^^-A= (M^-M + APj^i-M^-A, (87) 

where P is some penalty which should impose smoothness 
for b and c both angularly and as a function of time. For 
instan ce. For the b field we could use (see, e.g. IPichon et alJ 
<2002h l: 

P[b(m] = where Q(w) = (]S^„,]^), 

e J 

(so that large cv and t are less likely in the solution) and a 
similar expression for the c field which should also impose 
smoothing along velocities. The penalty coefficient. A, should 
be tuned so as to provide the appropriate level of smooth¬ 
ing. In practice, it might be necessary to impose further non¬ 
linear constrains on the solution, b, such as requiring that the 
excitation is locally as compact and connex as possible on the 
R 200 sphere. This can be done via some form of non-linear 
band pass filter in the prior, in order to limit the effective de¬ 
grees of freedom in b. 

Accounting for non-linearities. 

The non-linear solution, Eq. <5H may be formally rewrit¬ 
ten as a 2 = M 2 • b ® b, so that the perturbative inverse reads 

b 4 • A - • M 2 • • A) ® (M^^^^ • A), (88) 

(where M^ is defined by Eq. <871 1 provided the regime for 
the perturbative expansion applies. If not, we may still find 
the best non-linear solution to the penalized likelihood prob¬ 
lem of jointly minimizing ] ] A — M • b — M 2 • b ® b] ]^ -f AP, 
while using Eq. <88< as a starting point. 

When proper motions measurements are not available 
{i.e. we only have access to star counts), Eqs. <86l-i551 still 


5.4 Past history of galaxies : dynamical inversion 

Let us now see how the theoretical framework presented in 
Section |3 and Section may be applied to invert observed 


* Since our treatment of the dynamics (including the self -consistent 
gravity polarisation) is linear order by order, we may in principle 
recover the history of the excitation. 
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apply with some straightforward modifications, but the con¬ 
ditioning of fhe problem should decrease significantly, since 
the dynamics is less constrained. 

For a data set such as GAIA, we shall have access to the 
full 6-dimensional description of phase space for some of the 
stars (via radial velocity measurements and parallax) or at 
least 5 dimensional measurements (f, b, pif, fi},, Uy). 

Recall that in practice, the fields, b and c are respectively 
three-dimensional (2 angles and time) and 5 -dimensional 
(2 angles, time and 3 velocities). Consequently the in¬ 
verse problem is generically very ill-conditioned since data 
space is either two (£, h), four (f, }iP}, five -dimensional 

(£, b, fii, itb/Ti), or six dimensional {£, b, fig, pii,,7T,Vr). In fact it is 
anticipated that the conditioning is even poorer because the 
dynamical evolution involves damped modes, implying an 
exponential decay ( which corresponds to a major challenge 
for exfrapolation). It remains that the weakly damped modes 
should be tractable back in time up to some horizon, which 
will depend on the nature of the halo (via the conditioning of 
M;^ defined in Eq. |87|), the volume of data, and the signal to 
noise ratio in the measurements. 

Let us close this discussion of the inverse problem by 
emphasizing again the true complexity of the implementa¬ 
tion: Eq. <87> . and its non-linear counterpart, Eq. <88> . include 
via the dot product large sums over n and integrals over t. 
M and M 2 are functions of (which require a couple of 
integrals) and A^i, etc... which are themselves functions of K,- 
(Eq. <53H (which involves the underlying distribution func¬ 
tion, fo(I)/ srid fhe basis function, i/)["l(r) via Eqs iAfl . 142\ 
and j47|). 

Streamers (or tidal tails) in external galaxies may also 
be integrated backwards through the same procedure. It will 
involve the deprojection of the stream and of the underlying 
halo. 

In contrast to the Galactic inverse problem, it will in 
principle be possible to reproduce the inversion process on 
a statistical set of haloes, which would allow us to compare 
directly to the predicted statistical properties of the b and the 
c (though clearly the bias introduced by the penalised inver¬ 
sion would have to be accounted for). 

This completes our rapid survey of possible applications 
for the perturbative treatment of the dynamics of an open 
halo. 


6 CONCLUSION 

In the last few years, with the observational convergence to¬ 
wards the concordant cosmological model, a significant frac¬ 
tion of the interest has shifted towards smaller scales. In¬ 
deed it now becomes possible to project down to these scales 
some of the predictions of the model. This in turn offers the 
prospect of transposing there what has certainly been a key 
asset of modern cosmology, both observationally and theo¬ 
retically: statistics. This is a requirement both from the point 
of view of the (often understated) variety of objects falling 
onto an L* galaxy, but also because of the sheer size of the 
configuration space for infall. It is also a requirement from 
the point of view of the non-linear dynamics within the dark 
matter halo in order to account for the relative time ordering 
of accretion evenfs. was pointless to describe continuous in¬ 


fall, haloes are typically not in fully phase mixed equilibria, 
and the resulting fluctuation spectrum may seed or excite the 
observed properties of galaxies. 

In this paper, we aimed at constructing a self-consistent 
description of dynamical issues for dark matter haloes em¬ 
bedded in a moderately active cosmic environment. It re¬ 
lied entirely on the assumption that the s tatistics of the in- 
fall is w ell-characterized, as described in lAubert & Pichoi^ 
<2005albll . and that the mass of the infalling material (or to 
a lesser extend that of the fly-by) should be small compared 
to the mass of the halo. It also assumed that the halo was 
spherical and static or evolving adiabatically (Section|4}. The 
emphasis was on the theoretical framework, rather than the 
details of the actual implementation. In other words, we 
aimed at describing a self-consistent setting which allows us 
to propagate the cosmological environment into the core of 
galactic haloes. 

In Section 13 we derived the dynamical equations gov¬ 
erning the linear evolution of the induced perturbation by di¬ 
rect infall or tidal excitation of a spherically symmetric (inte- 
grable) stationary dark matter halo. The simplified geometry 
of the initial state allowed us to focus on the specificities of an 
open system. Specifically, we revisited the influence of the ex¬ 
ternal perturbations on the spherical halo, and extended the 
results of the literature by considering an advection term in 
the Boltzmarm equation. This approach was compared to the 
classical Green solution in Appendix^ Note that both the 
intrinsic properties of the halo, via the distribution function, 
F (Eq. I16H. and the environment, via (s*^, ip^) (Eq. <17H. of 
the infall and the tidal distortion were accounted for. Clearly 
the subclass of problems corresponding to tidal perturba¬ 
tions only will turn out to be easier to implement at first. Ap- 
pendix^presents the details of the angle action variables on 
the sphere together with an explicit expression for the kernel, 
K, and carried out a test case implementation of the statisti¬ 
cal propagation of an ensemble of radial excitations with a 
powerspectrum scaling like 

In Section 13 we derived the non-linear response of the 
galactic halo to second-order (Eq. <44H in the perturbation 
(and to order n in Appendixloltogether with the correspond¬ 
ing N-point correlation function) to account for tidal strip¬ 
ping and dynamical friction. The dynamics was "solved" it¬ 
eratively, in the spiri t of th e successful approach initiated 
in cosmology bviFrvI il984h and considerably extended by 
bernardeaul il 992ft . In particular we presented and illustrated 
a set of diagrams (Fig. |3, lUlH . each corresponding to the 
contribution of the perturbation expansion. Though the ac¬ 
tual implementation of the non-linear theory is going to 
be CPU intensive, we argue that it will improve our im- 
derstanding of the competing dynamical processes within 
a galactic halo. In particular, we discussed how this explicit 
theory of non-linear dynamics provides the setting in which 
substructure evolution (and destruction) will have to be car¬ 
ried, in order to account for e.g. tidal stripping. 

In Section 21 we presented the Fokker-PIanck equation 
governing the quasi-Iinear evolution of the mean profile of 
the ensemble averaged halo embedded in its cosmic environ¬ 
ment. Specifically we showed how the infall, driff and dif¬ 
fusion coefficients (Eqs. <G13l - td^ ') are related to the two- 
point correlation of the tidal field, and incoming fluxes. Ap- 
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pendix|3gives a derivation of this equation from first princi¬ 
ples, while in the main text, we focused on the bibliographic 
context and possible applications. The key physical ingre¬ 
dient behind this secular evolution theory was the stochas¬ 
tic fluctuation caused by the incoming cosmic substructures. 
The key technical assumption was that the two time scales 
corresponding to the relaxation processes and the dynami¬ 
cal evolution decouple. Hence we could assume a hierarchy 
in time so that the distribution function is constant in time 
when computing the polarization. 

Finally, in Section|^we considered in turn a few classi¬ 
cal probes of the large-scale structures which had been used 
in the past to constrain the main cosmological parameters 
and the initial power spectrum, which we transposed to 
the galactocentric context. Note that these are built upon ob¬ 
servables, hence they may be used to constrain the boimd- 
ary power spectrum of the fln- Since Eq. I65> , I71t . <EH . 
and mn) involve different combinations of (flit^n')' they 
will consfrain them at different scales with different biases, 
which should ultimately allow us to better characterize the 
power spectrum. This situation is the direct analog of fhe cos¬ 
mic situation, where the different tracers (weak lensing, Ly-a 
forest, CMB etc..), constrain different scales of the cosmolog¬ 
ical power spectrum (with different biases). Note also that 
our knowledge of the statistical properties of the boundary 
(via the bn and Cn coefficients) together with some assump¬ 
tions on the equilibrium Fg all ows us to generate given real¬ 
izations of fhe fln as shown in lAubert & Pichc^ ^200.5a^ and 
therefore virtual observables for any of fhese data sets, for 
the purpose of e.g. validating inverse methods. We investi¬ 
gated the consequences of the infall down to galactic scales 
and showed how it could be used to account for the observed 
distribution of disk properties (spiral winding, warps etc... 
). We demonstrated how the analytical model (both linear 
and non-linear) are quite useful when attempting to'Tnvert" 
the observations for the past accretion history a given galaxy. 
This stems mainly from the fact that perturbation theory pro¬ 
vides an explicit scheme for the response of the system, in 
contrast to the algorithmic procedure corresponding to N- 
body simulations. 

Again let us emphasize that Eq. iTsl and its non-linear 
generalization ID3I and IC7I yield in principle the detailed 
knowledge of fhe full perturbed distribution (inside f?20o) 
later times. (This is to be contrasted to the situation in N- 
body simulations where the response of the system is par¬ 
tially hidden by the mean profi le of the halo , which requires 
first identifying substructures iAubert et alJ i20n4l ')'). Hence 
we should be in a position to weigh the relative importance 
of the environment (via and ip'^) against the inner proper¬ 
ties of the galaxy: the unperturbed distribution function of 
the halo, F(I), (its level of anisotropy, the presence of a cen¬ 
tral cusp etc... ) the disk, (its mass, its profile, its distribution 
function, F(J) etc... ). 

The work presented here derives from the fact that it 
was rea l ized that the biorthogonal projection pioneered by 
iKalnaij Jl97(j) could be applied order by order to the per¬ 
turbative expansion of the dynamical equations. Yet this in 
turn required the knowledge of the relative phases involved 
in the perturbation, which involves characterizing the prop¬ 
erties of the perturber. The characterization only made sense 


st atistically in o rder to retain the generality of the approach 
of iKalnaia ll97(h . Hence the emphasis on statistics. 


6.1 discussion & Prospects 

Our purpose in this paper was to address in a statisti¬ 
cally representative marmer dynamical issues on galactic 
scales. We also advocated using perturbation theory in angle- 
actions in order to explicitly propagate this cosmic boundary 
inwards in phase space. As was demonstrated in the paper 
(and shown quantitatively in Section IbH . this task remains 
in many respects quite challenging. 

One of the limitations of the above method is the re¬ 
liance on numerous expansions combined to the special care 
required in their implementation. One could argue that this 
level of sophistication might not be Justified in the light of fhe 
weakness of some of the assumptions. Indeed, we are limited 
to systems with spherical geometry whereas galaxies most 
likely come in a variety of shapes. This assumption could be 
lifted provided we compute the modified actions of the flat¬ 
tened spheroidal equilibrium using perturbation theor y for 
the equilibrium in the spirit of iBinney & Spere^ il984h . but 
implies a higher leyel of complexity; (it would also require 
statistically specifying the orien tation of the halo re latiye to 
the infall, as discussed in part in lAubert et alJ 12004) 1. We as¬ 
sumed here that the perturbation was relatiyely light, which 
excludes a fraction of cosmic eyent which might dominate 
the distribution of some of the obseryables. 

Section lEl Section |5l and Appendix El presented a few 
po ssible application s for the framework described here and 
in lAubert fe Pirhorj l2005all . These galactic probes would 
need to be further inyestigated, in particular in terms of ob- 
seryational and instrumental constraints. The biasing spe¬ 
cific to each tracer should be accounted for. The second-order 
perturbation theory needs to be implemented in practice to¬ 
gether with the diffusion coefficients of Sectional following 
Section |B| and extending Section IbII Similarly, the identifi¬ 
cation and eyolution of substructures within t he halo men¬ 
tioned in Section l3.3.1l deseryes more work. In lAubert et alJ 
l2nn4l . we showed that the accretion onto L* haloes was 
anisotropic; the dynamical implication of this anisotropy will 
require some specific work in the future. 

We will need to demonstrate against N-body simula¬ 
tions the releyance of perturbation theory for dynamical fric¬ 
tion; in particular, we should explore the regime in which the 
second-order truncation is appropriate, and at what cost ? 
Note that truncated perturbation theory implies that modes 
will ring foreyer. At some stage, one will therefore haye to 
address the problem of energy dissipation. 

Implementing a realistic treatment of the infalling gas 
will certainly be amongst the more serious challenges ahead 
of us. This is a requirement both from the point of yiew of the 
dynamics but also from the point of yiew of conyerting the 
aboye predictions into baryon-dependent obseryables. The 
description of the gas will require a proper treatment of the 
yarious cooling processes, which can be quite important on 
galactic scales. In particular, the thickening of galactic disks 
is most likely the result of a fine-tuning between destructiye 
processes such as the tidal disruption of compact substruc¬ 
tures on the one hand, and the adiabatic coplanar infall of 
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cold gas within the disk. In fact, the nonlinear theory pre¬ 
sented in Section 13 and Section could be extended to the 
geometry of disks to account for the adiabatic polarization 
towards the plane of the disk. 

Note that we assumed here that transients correspond¬ 
ing to the initial conditions where damped out so that the 
response of the system was directly proportional to the exci¬ 
tation. The underlying picture is that of a calmer past, which 
in fact is very much in contradiction with both our measure¬ 
ments and common knowledge on the more violent past ac¬ 
cretion history of galaxies. Indeed, infalling subclumps will 
contribute via the external tidal potential at some earlier 
time, and the larger the lookback time, the relatively stronger 
the importance of the perturbation (since the intensity of in¬ 
fall is in fact an increasing function of lookback time). We 
are therefore facing a partially divergent boundary condi¬ 
tion. Because of the characteristics of hierarchical clustering, 
the actual bootstrapping of the analytical framework is there¬ 
fore challenging. This could be a problem in particular for 
non linear dynamics, where the coupling of transients may 
turn out to be as important as the driven response. The im¬ 
portance of these shortcomings will need to be addressed in 
the future. 

Finally let us note that the theory described in Section|3 
Section |3 and Section |4] describe perturbative solutions to 
the collisionless Boltzmann Poisson equation in angle-action 
variables, and as such are not specific to the description of 
dark matter haloes. It could straightforwardly be transposed 
to other situations or geometries provided the system re¬ 
mains integrable. As mentioned in Section lZSl the stellar dy¬ 
namics aroimd a massive black hole would seem to be an 
obvious context in which this theory could be applied. For 
instance, we might want to investigate the capture of streams 
of stars by an infalling black hole. In a slightly different con¬ 
text, note in passing that the above theory could also be ap¬ 
plied to celestial mechanics, since an angle-action expansion 
corresponds to an all eccentricity scheme. 

Let us close this paper by a summary of the pros and 
cons of the theory presented here. 

Possible assets: 

• fixed boundary: localized statistics; 

• fluid description: no a priori assumption on the possibly 
time- dependent nature of the objects; 

• non-linear explicit treatment of the dynamics: proper ac¬ 
count of the self-gravity of incoming objects and statistical 
accounting of causality; 

• dynamically-consistent statistically-representative treat¬ 
ment of the cosmic environment; 

• customized description of resonant processes within the 
halo via angle action variables of universal profile; 

• ability to construct one- and two-point statistics for a 
wide range of galacfic observables. 

• theoretical framework for dynamical inversion and sec¬ 
ular evolution 

Possible drawbacks: 

• weak perturbation w.r.t. spherical stationary equilib¬ 
rium: not representative of e.g. equal mass mergers; 

• complex time dependent 5D boimdary condition; 


• ad hoc position of the boundary; 

• no obvious truncation of two-entry perturbation theory; 

• no account of baryonic processes; 

• inconsistency in relative strength of merging events ver¬ 
sus time; 

• non-Gaussian environment probably untractable; 

• finite temporal horizon given finite imax) 

• no statistical accounting of linear instabilities. 
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Figure 5. diagrammatic representation of the expansion to second- 
order given in Eqs. 14.31 - 1441 . The top diagram states that one should 
sum over all orders in the coupling in order to model the non-linear 
response of the halo; The second diagram from the top stands for 
Eq. 1431 and the third for Eq. 1441 . The loops correspond to the self¬ 
coupling i.e. the self-gravity response of the halo to the perturbing 
flow. The second diagram corresponds to the "propagation" of the 
double excitation (see also AppendixIXlfor a discussion of the distri¬ 
bution function propagator in angle-action variables): the input are 
the external potential, >p‘ (through its bn coefficients) and the source, 
s‘ (expanded over the Cn coefficients); the output is the coefficient of 
the expansion of the inner potential. The coupling is achieved via the 
operator fC; and Q, defined by Eqs. 1421 - 13.l.Tl and ID4I . while the 
contraction is achieved by Eq. 1451 . 1461 and is represented by the 
wiggly horizontal line. 
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APPENDIX A: LINEAR PROPAGATOR IN ACTION ANGLE 


As mentioned in Section l2.4.1l it is useful to regard the open collisionless system as a segmentation of the source for the prop¬ 
agator (ie the Green function of the coupled Poisson Boltzmann equation), where one distinguishes two contributions for the 
initial distribution: the contribution at R 200 with Vr < 0 (what we describe as the source term in Eq. and the contribution 
beyond _R 200 or at _R 200 with Vr > 0 (what we describe as the tidal field in the main text). In order to make this comparison, let us 
derive generally (without any reference to a boundary for now) the Green function satisfying the linearized Boltzmann Poisson 
equation. Let us call G(w, I, t|w', I', t') this Green function; it obeys: 


9G „ ^ „ 

-l- cu • VwG + • Vv 

dt dl 


/ 


dr" 


j-dv"G(r", v", t|w',l', t') = <5d(w — w')<5d(I — l')i5o(t — t'), 


(Al) 


so that the distribution fimction at (I, w, t) reads 


/(w, I, t) = J j J dI^G(w, I, t|w', I', t')/(w',E, f'). 

Let us define the linear propagator, iiai,k,wo(I|lO/ 


^ (k.cu~a>) ((1 - )„, n ' 5 (k'.cu^-a;) ^ 






(A2) 


(A3) 


so that the distribution function, /(L w, t), at time t, with action I and angles w induced by the propagation of the distribution 
at earlier time t', with action I', and angles w' reads 

/(I,w,t) = J dt' J dl' J dv/'Y^ y da;exp(za;[t-t'] - !k'[w-w'])!i^^k,wo(I|I')/(P.w',L). (A4) 

It is interesting to contrast Eq. US to the propagator found bv llchimarul il973h for the uniform plasma. In particular, the gradient 
of the density profile breaks the stationarity in w — wq of the propagator, Eq. <A3> . Note also that the first term on the r.h.s. of 
Eq. <A3I corresponds to free streaming inside the halo (i.e. dark matter particles describing their unperturbed orbits), and reads 
in real space 

Gfree(I/w, t|l', w', t') = YY [ dmexp(z6t;[t — t'] — !k • [w — —LI = < 5^(1 — I'),5q(w — W — w[t — t']), 

^ J k • cu — cu 

while the second term in Eq. <A4I corresponds to the self-gravitating polarization of the halo induced by the perturbation. 
Note that since the field dynamical equation is solved with a right-hand side {i.e. a source breaking the mass conservation 
in phase space), Liouville's theorem is not obeyed anymore: a new fluid is injected into the halo. We may now assume that 
in Eq. <A2> . = /(w',I', t') is split in two: one contribution from dark matter particles exiting K 200 or beyond R 200 } 

another contribution describing particles on R 200 with negative radial velocity. The former component may then be resumed 
over the corresponding region of phase space with a 1 /1 r — r' | weight, and yields ip". The latter corresponds to (t'). 


APPENDIX B: IMPLEMENTATION 

Let us describe in this appendix in greater details how Section|2lare implemented in practice, while focusing here on a simple 
isotropic halo {i.e. F(I) = F{E), where E = /I 4 - Yjr) is the energy, and Y(r) the unperturbed potential). We will show here 
how to comp ute the operator, K, (defined bv Eq. <16H and elements of O (defined by Eq. II71 L for the corresponding basis, 
following e.g. iTremaine & Weinberel il984h . iMuralil h999ft . ISeguin & DuprazI il994h . We will then implement in practice the 
average induced correlation triggered by some ad hoc colored radial perturbation. Similarly, one could compute the non-linear 
coefficients, |[ ] entering Eq. 1.391 . but the implementation of the non linear formalism of Section|3and Section01are beyond the 
scope of this paper. 


B1 Detailed angle-action linear response for isotropic spheres 

The three-dimensional nature of galactic halo makes the implementation slightly more complicated than one would think at 
first sight. The assumption that the halo is spherical allow us to assume that the equilibrium is integrable. Hence the action 
space is effectively at most two-dimensional, but configuration space remains three-dimensional (though one angle is mute). In 
practice, this implies that integration over action space, occurring in e.g. Eq. 1411 is effectively two-dimensional. On the other 
hand, the sum over k involves three indices, each corresponding to a degree of freedom. 

Let us define p as the radial action, l 2 = L as the total angular momentum and I 3 = as the z-component of the angular 
momentum, so that 

h = - f‘ drJ2[E~'¥{r)]-ll/P. 

^ Jr, 
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Here and rp are respectively the apoapses and periapses of dark matter particles. This defines I = (fi, I 2 , 13) introduced in 
Section lZn Similarly, Let us define the corresponding angles, w = (zvi, ZV 2 , M’s) (see Fig. 0) given by: 

dr(a; 2 -f 2 /''^) 


ZVi = OJi 


f 

Jr„\ 


dr 


V2[E-Y(r)l. 
where cos ^ = Lz/L. 


rr{\.,wi) 

==, W2{l,Wl)=X- —r 

f2/r2 Jryi) y; 


2[£-¥(r)]-f2/,2 


W 3 = (p — asin(cot(/3) cot(0)), (Bl) 


B2 Computing the linear response operator 

Following very closely the notation nf iMur inimi, let us introduce a bi-orthogonal basis constructed around spherical har¬ 


monics: 


| 0 (hO= I] and Xp{t,t)= ^ , 

£mn £mn 

for respectively fhe densify and the potential. IWeinbe^ il989ll suggests the following potential-density pair 
R^^^^ie{oinr/R), and df”'(r) = -p 


//(r) = - 


47tG\/2 

®n|yr(fl^n)| 


n\/2 

lieMl 


R-^^^idc^nr/R), 


(B2) 


(B3) 


where ji stands fo r the spherical Bessel functio n and where obeys the relation an/f-i/ 2 (“n) = 0. Here R is the truncation 
radius of the basis. IS;rnquist & Ostrike'3il992h suggest another set of (non-normalised) biorthogonal functions defined by : 


uf (r) = - 


(1 + r)2Hl 


and C(r) = ^ 


J-i 


2 n (1 + r)2H3 




where K„( = nl2{n + 4f + 3) + (f -I- l)(2f + 1), f = (r — l)/(r + 1) and C^(x) stand for ultraspherical polynomials. 
The action angle transform of the potential basis is given by : 


Wl"{l) = xp^{l) = ^ J dwiexp{-ikiWi)ufdr)exp[tk2{x-W2)], 

We may now rewrite Eq. <16> as 


K"'(T-f) = -4C' 


TttG 


// 


dE-^^ wexpfik-a;(T- t)]W^^”'(I) 

1 k 


Km + 


471 




where 

pi- = J drrHl"‘{r) 


dr ' 


and where 

2^k2-l (^l_k2)[Yfl/2{e + k2 + l)] 
7t 2 (^-Hfc2)!F2[l/2(f-fc2) + l] ' 


if £ + k 2 even, else 0 . 


(B4) 


(B5) 


(B 6 ) 


(B7) 


(B 8 ) 


Here T is fhe sfandard Gamma fimction. Note that Xk(I) accounts for the fact that the response is computed in a non inertial 
referential frame. To take into account the barycentric drift of fhe halo, the perturbed Hamiltonian should include the induced 
inertial potential • r, where a;, is the acceleration of fhe barycenter in the frame of the unperturbed halo. Its action-angle 
transform is given by : 


1 f" 

Xk(I) = ^ J dwiexp{-tkiwi)rexp[tk2{xp ~W2)]. 


(B9) 


As can be seen from Eq. <B 6 > . this inertial contribution is limited to the dipole component (f = 1) of the response : as expected, 
it is equivalent to a spatially homogeneous field force.® 


® Technically speaking, the Sn dependence arise from the fact that r is expressed as a function of Yim (O) spherical harmonics. 
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B3 Implementation and validation 

The actual computation of the linear response of the halo 
to a tidal field is a two-steps procedure. First the kernel K 
must be computed via Eq. <B6I . It involves an integration 
over the orbits' space and requires to Fourier transforms the 
the biorthogonal basis (W and X quantities) along orbits. It 
can be done by "throwing orbits" in the equilibrium potential 
and finding the associated sets of (I, w) in the halo's model: 
such a procedure provides the angle dependance of the basis' 
functions for a given action. Knowing W(I),X(I), cu(I) over 
a given sampling of the I space, Eq. <B6I can be computed. 
In order to achieve high computing efficiency and accurate 
responses, we implemented the calculation of Eq. IB6> in a 
parallel fashion, where the integrals in each-subspace of the 
action space are computed by a different processor. 

Second, the expansion a(i) of the halo's response is com¬ 
puted eithe r by iteration or b y means of a Volterra's Equation 
solver (e.g. lPress et alJ il992h ~). We found that both methods 
give very similar results and differ only by their time con¬ 
sumption. The iterative method can be very fast if a proper 
initial guess is available but if it is not the case it may take 
a significant amount of time to achieve convergence. Con¬ 
versely, the Volterra solver's time consumption is fixed for a 
given time resolution. 

In order to validate our impl ementation, we s et up two 
tests. The first one is suggested bv lWeinberd il989h . A Plum¬ 
mer's halo is embedded in an homogeneous force field and 
should experience a global drift decribed by the potential's 
response: 

!/;(r,f) = -rb(f). VOW, (BIO) 


where stands for the barycenter position and O stands for 
the equilibrium potential. We chose the force field to have 
a flosin(vf) time dependance with ag = 0.01 and v = 0.01. 
The Plummer model has a unit mass M and caracteristic ra¬ 
dius b. The response was computed using a 60 x 60 sampling 
in (£, L) and 20 radial terms of the basis given by Eq. <B3t . 
We switched off the drift compensation modelised by the X 
term in Eq. <B6t . Fig. <BH shows the response computed at 
t = 10 (in units of yW/GM) along with the prediction given 
by Eq. IBlOt . Clearly the two responses coincide, providing a 
first validation of our implementation. 

A second test involves reproducing the contraction of a 
Hernquist's halo induced by a central spherical mass (which 
would model the presence of a galaxy for example). This cen¬ 
tral mass is assumed to follow a Hernquist's profile, whose 
potential is given by: 


®w 


GM 
r + a' 


(Bll) 


The halo has a unit mass M and caracteristic radius a, while 
the central object has a final mass of nip = 0.001 and 
flp = 0.25 as a constant caracteristic radius. The perturber is 
turned on at f = 0 and follows a mp{t) = mp(t/1 jp — 
2 {t/tj:p) temporal evolution, where iy is the final time step. 
We compare the linear response at f = iy with the simula¬ 
tion of the same test-case using a perturbative particle code 
(Magorrian, private communication). The response was com¬ 
puted using a 60 X 60 sampling in 13 subregions of the whole 
(E,L) space and 21 radial terms of the basis given by Eq. <B4> . 



Figure Bl. Isocontours in the X-Y plane of the potential's response of 
a Plummer Sphere embedded in a homogeneous force field (see text 
for details). The force field is aligned along the X axis. The dashed 
line stand for the prediction and the solid line stands for the linear 
calculation presented in this paper. 


Self gravity of the response is not taken in account in both 
methods. Fig. IB2I shows the comparisons between the two 
type of calculations, made for two different growing time fy. 
Clearly the two methods are in good agreement. One can see 
that matter is dragged toward the center and the longer it 
takes to grow the perturber, the further from the center are 
the affected regions. 


B4 Statistical propagation: a test case 

In this section, we compute the two-point statistics of a halo 
responding to a simple type of tidal perturbation as an illus¬ 
tration of statistical propagation. Without any assumption on 
the type of perturbation, we recall that the two-point statis¬ 
tics of the halo's response is given by Eq. <29> and can be 
derived directly from the perturbations' statistics. Let us sim¬ 
plify the correlation's computation by assuming that the halo 
is only tidally perturbed, so that Eq. 1291 reduces to : 

(^a • = ( [k • b] • (1 - K)-i • (1 - K)-!*^ • [k ■ b] ).(B12) 

Furthermore, let us also (rather crudely) assume that the tidal 
field is monopolar, and has a radial dependence equals to the 
Nth element of the radial basis which diagonalize the Pois¬ 
son equation. Then, the tidal perturber's coefficient can be 
written as : 

= b{t)SnN^t0^m0, (B13) 

where the perturbing tidal field is described by : 

f{r,at) = b{t)u^g{r). (B14) 

Since no radial coupling occurs, the halo's response can be 
simply written as : 
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Figure B2. The radial profile of the density response of a Hernquist's 
halo due to a central pertuber. Lines stand for the linear response 
of the halo as the central pertubrer grows over a t = Q.63t^ (dotted 
line) and t = lOfd timescale (dashed line), fd is the dynamical time 
of the main halo within its core radius. Radii are given in units of the 
main halo's core radius. Density units are in code units. Superim¬ 
posed are the calculations of the same response using a perturbative 
particle simulation (Magorrian, private communication). No scaling 
has been applied and the two methods agree quantitatively. 

xp{r,d,t) = a{t)u^Q{r). (B15) 

Consequently, the only remaining degree of freedom is the 
temporal variation of the tidal field. If we consider an ensem¬ 
ble of tidal environments and if we assume stationarity and 
gaussianity of the induced perturbations, it will be described 
by the temporal two-point correlation function of hit) coef¬ 
ficients, or equivalently by their temporal power spectrum 

PbiH 

Pi,{v) = {HDb*{v)), (B16) 

where v stands for the frequency. If the temporal power spec¬ 
trum of the response is given by Pa (v) then Eq. <B12> reduces 
to: 

P„(t/) = («(!/)«* (v)) = -^%f^Pb(v). (B17) 

Eq. IET7\ simply states that the frequency structure of the 
'tidal noise' is transmitted to haloes via a (scalar) transfer 
function given by the response kernel. 

Let us furt her describe our test halo again by a Hern¬ 
quist's model iHernquisJ il99ril) ). The corresponding ker¬ 
nel is computed fol lowing the procedure describ ed by Sec¬ 
tion IB2I using the iHernquist & OstrikeI^ h 992ft potential- 
de nsitv pair (see Eig. Eurthe r deta ils can be found 

in iMuralil il 999ft . ISeguin & Dupra j il994ft . The radial de¬ 
pendence of the tidal perturber is given by the N-th po¬ 
tential function uw(r) = u9?(r) of the basis described by 
IHernquist & OstrikeI^^1992ft . The associated density function 
is givenby d]y(r) = = Aw]y(r)/47rGandexamplesof 

such profiles are given in Fig. <B5I along with the halo's pro¬ 


aoi 0.1 1.0 10.0 100.0 0.1 i,o io.o looo 



Figure B3. The halo's density profile chosen for the statistical prop¬ 
agation's example follows a hernquist model (top-left panel). We 
apply a monopolar tidal field >p‘^(r) with a radial structure given 
by the d^{r) function of the Hernquist and Ostriker biorthogonal 
basis. Here are shown the corresponding density profiles p‘ir) = 
AtpW4:nGforN = 1,5,10. 

file. For simplicity, the tidal frequency distribution has been 
chosen to follow a power law: 

Pfc(v) ~ v-^. (B18) 

This power law describes the ensemble frequency behavior 
and therefore a single realisation of the tidal noise may de¬ 
viate from this relation as long as statistical convergence is 
achieved. Fig. <B4> shows both an example of the time de¬ 
pendence of such a perturber and the time dependence of the 
induced response. One can see that the halo acts a low pass 
frequency filter and do not recover all the high frequency fea¬ 
tures present in the tidal field. Also, the halo response ap¬ 
pears as delayed in time, reflecting the effect of the halo's 
own inertia. 

The same computation was performed for an ensem¬ 
ble of 1000 different tidal perturbations. Fig. <B3I shows 
the power spectrum Pi,(v) averaged over all the realisations 
along with Pa{v) averaged over the 1000 haloes' responses 
(shown as symbols with error bars). Pi,(v) departs from a 
power law at low frequencies (v < 50 in code units) because 
of the finite time range over which the tidal field is applied 
(not shown here). At higher frequencies, the perturbers' fre¬ 
quency distribution follows exactly Eq. IB18> . Independantly, 
Pfl(v) is directly predicted from Pi,{v) using Eq. iB17> . with¬ 
out relying on the computations of individual responses, and 
shown on the same plot as solid lines. Clearly the predicted 
power spectrum of the response matches the statistically av¬ 
eraged one and even reproduces 'bumpy' features seen at 
various frequencies. The filtering effect of the halo response 
can still be seen in the predicted spectra : Pa{v) follows the 

law at low frequencies but exhibits a steeper slope at 
higher frequencies. This cut-off effect is more important for 
large scale perturbations (low N) and reflects the fact that 
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Figure B4. An example of time evolution of the tidal field's ampli¬ 
tude b{t) (plain line). Its power spectrum Pi,{v) follows a law. 
Assuming an N=2 radial dependance, the amplitude fl(f) of the in¬ 
duced halo's response can be computed (dashed line). The halo does 
not respond to high frequency features and globally its response is 
slightly delayed, reflecting its own inertia. 

perturbations at high frequencies are unable to 'resonnate' ef¬ 
ficiently with the halo's large scale modes. Conversely, tidal 
perturbations with features on small spatial scales (large N) 
are more likely to induce large frequencies and preserve 
the frequency structure of the perturbation. Moreover, the 
'bumps' seen in the Pa{v) curves reflect the eigen frequencies 
of the halo. The pertuber's scale-free spectrum hits reson- 
nances which react in a stronger fashion than any other fre¬ 
quency. Again, these resonances occur at larger frequencies 
as the radial order N increases : smaller radial scale pertur¬ 
bations relate to shorter caracteristic time scales. 

The illustration presented in this section is admittedly 
simplistic but hints at the possibilities which can be foreseen 
for statistical propagation : for a given set of constrained en¬ 
vironment, predictions on the statistics of the induced re¬ 
sponse can be made without relying on the computation 
of individual realisations. Predictions on spatial or spatio- 
temporal correlations of the halo's response can be made fol¬ 
lowing the same procedure. It will possibly allow us to study 
the impact of the different scales of accretion or potential, 
the influence of the rate of change of these pertubrations and 
their relative relevance on the statistical properties of matter 
within the halo as discussed in Section|3 



frequency v 


Figure B5. An example of statistical propagation. The average power 
spectrum of the 1000 tidal perturbations applied to the Hernquist 
halo (top curve) follows a law (dashed thick curve). Symbols 
stand for power spectrum of the halo's response averaged over 
the 1000 realisations of tp‘{r,t) with four different radial depen¬ 
dences (with N = 1,3,5,10, from top to bottom). The superimposed 
curves show the direct predictions on the power spectra, following 
Eq. IR17I . For clarity these curves have been divided respectively by 
1, 20, 40 and 70. Frequencies are in code units. 
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APPENDIX C: SECULAR EVOLUTION WITH INEALL 

Let us derive in this appendix the secular equation for the evolution of the ens emble average hal o em bedded in a typical cosmi c 
environment (with infall and tidal field). This follows the pioneering work of lWeinberd ^2001a^ and iMa & BertschmgerH2004ll. 
The settings in which they derive their coefficients differ: their starting point is the kinetic closure relation given bv l^^montovichl 
h967tl . llchim^ il973h . lGilberl il97nll who note that the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy may be 
closed while assuming that the two-point correlation function will relax on a shorter dynamical time scale, whereas the one point 
distribution function evolves on a longer secular time scale^^. Hence if one assumes that the distribution function F entering 
the linearized equation, Eq. can be considered to be constant, then the second-order equation in the BBGKY hierarchy is 
automatically satisfied while the r.h.s. of fhe first equation is proportional to the propagated (via Eq. <A3H excess correlation 
induced by the dressed clumps. This kinetic theory has been successfully applied in pla sma physics, lea ding to the so called 
Lenard-Balescu ihenardl h96lll . lB'alescul h96.'jB collision term, and was also transposed bv lWeinberd il99.dl for a multi-periodic 
"stellar" system. 

Note that the BBGKY hierarchy is a 1/N expansion, where N is the number of particles in the system. Formally it would 
make sense here to identify N as a measure of fhe dumpiness of fhe medium, buf fhis definifion is qualitative only. We rely 
here on the same time ordering hierarchy, but the degree of dumpiness in the system is explicitly imposed by the boundary 
condition. In this appendix, the derivation is carried from first principles, while relying on an explicit infall and tidal field. 


Cl Quasi linear equations in angle action variables 

The collisionless Boltzmaim equation of an open system may be written as: 

7,2 

+ {H,F} = S" + s", with H = — + (Cl) 

at 2 

where F is defined by 

F(I,w,t,T) = F(I,T) 4-/(1,w,t), and Y(I, w, t, T) = To(I,w, T)+!/;(!, w, t)-f i/>''(I,w, t), (C2) 


wifh F describing the secular evolution of the DF and / describing the fluctuations of the DF over this secular evolution. In 
Eq. ICII . the r.h.s. stands for the incoming infall, both fluctuating (s'^[I, w, f]) and secular (5*^(1, T)). Since this system evolves 
secularly because of ifs environment, these actions are not conserved. The last two term on the r.h.s. of Eq. tC.2\ represents the 
fluctuating component tracing the motions of clumps within the environment of fhe halo^^. Note that since F(I, T) is assumed 
to depend here only on the action, it represents a coarse-grained distribution function (averaged over the angles) for which we 
make no attempt to specify where each star is along its orbit nor how oriented the orbit is. Note also that the canonical variables 
I and w are the actions and the angles of the initial system. Developing the collisionless Boltzmaim equation, Eq. O. over the 
secular and the fluctuating expansion leads to: 


dt dt dw dw ai au 


dF 

3w 31 


3U 31 31 3w 


S" + sL 


(C3) 


This equation involves two time scales, t and T. On the fluctuation time scale, t, secular quantities can be described as static, 
leaving only the linearised open collisionless Boltzmann equation (Eq. j^): 


dt 3w ^ 3w 3w 31 


(C4) 


where the amplitude of / is of firsf-order compared to F and involves only the fluctuating part of the external forcing, (I, w, t). 
On a longer time scale, T, the Boltzmann equation, Eq. O can be T-averaged, considering that the average of fluctuations are 
zero on such time scales. This leads to a second equation: 


m 

3T 


'dip dipn 3/\ 

3w 3w 31 


dip dip‘^ 

ST “Sf 


3w 




(C5) 


The brackets denotes averaging over a time longer than the typical time scale of fluctuations: 

l-T+AT/l 

{Y)t = 1/AT / dtY{t). 

Jt-AT/2 

The time interval, AT, should be chosen so that a given dark matter particle describing its orbit will encounter a few times the 
incoming clump at various phases along its orbit. Because the incoming clump is subject to dynamical friction, the resonance 
will only last so long, and induce a finite but small kick, AI during AT. Because the infall displays some degree of temporal 
and spatia l coher e nce, we may not assume t hat the successive kicks are uncorrelated, in contrast to the situation presented by 
IWeinberd <2001ah . lMa & Bertschinged i2004l or the classical image described in Brownian motion. In other words, when we 


This time ordering is originally due to lBogolvubov & GurovI Jl947ll 

We neglect here the secular drift of the external potential which should slowly shift the frequencies, u> in Eq. 01 
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write an effective microscopic Langevin counterpart to the corresponding Fokker Planck equation, it will involve a coloured 3D 
random variable (see Eq. <601 1. 

Derivative and averaging may be exchanged considering that F and Y evolve slowly with respect to time. Terms involving 
the product of two first-order quantities survive to the time averaging because we cannot presume that the response in distribu¬ 
tion function and potential within K 200 uncorrelated. In order to evaluate those quadratic terms, we may integrate Eq. <C14> . 
while assuming that F(I, T) is effectively constant w.r.t. time t. The solution, Eq. J7J, may then be reinjected into the quadratic 
terms in Eq. E3 so that they involves terms such as 

It • 1^ = - ( E ® k2j : ^ - 

E (I, f)k2 • ^ f (k (C6) 

ki,k2 J CO j 

where we may factor the action derivative of F(I, T) out of the r-time integral because the secular distribution is assumed to 
be constant over a few dynamical times. Since the l.h.s. of Eq. IC5I does not depend on w, we may avera ge its r.h.s. over dw. 
This i mplies that in Eq. EH, only the kj = — k 2 terms remain. We rely effectively on the averaging theorem isinnev & Tremain3 
to convert orbit averages into angle average. The corresponding evolution equation hence depends on the actions only, 
as expected. Note that in doing so, we assume that no other resonances matter. The secular equation, Eq. <C5< . becomes finally 
after some similar algebra for the other contributions:^^ 

^ u ^ c ^2 u 

- = (Do(I)) - (Di(I)) . — - (D 2 (I)) : (C7) 

where 

(Do(I)) = ^ j (S,)rdw+^^;k. E t) + xpi* {1, t)] e*'‘^("-')s,(k,I, rjdr) ^ , (C8) 

while the drift coefficient, Dj, obeys 

(Di(I)) = ^£kk.E (^[V^i;(I,f)+C(kf)] + 

and the diffusion coefficient, D 2 , is given by 

(D2(I)) = ^^k®k[V;^(I,t) + C(kf)] + . 

Note that the infall coefficient, Dg, includes both the secular infall, and a contribution arising from the possible correlation 
between the fluctuating tidal field and the fluctuating infall. It may be an explicit function of time, T, reflecting the fact that, 
as more mass is accreted, the profile of dark matter changes with time. The coefficients Dg, Dj and D 2 are also an implicit 
function of time because of the time average, ( )t and via the secular distribution function, F(I, T) which occurs in t) 
through Eq. <18< . Clearly, if the potential, and/or the source term are completely decorelated in time, so that (!/’k(I, t)ip^{l,T))j oc 
^D(f - t) and t)s^*(I,T))T “ t), Eq. <C10t or would vanish. Provided AT is long compared to the typical 

correlation time of the potential (and/or the source term), we may take the limit t —» 00 in the integrals entering Eqs. <C8< - <C10> . 
Note finally that Eq. <C7< does not derive from a kinetic theory in the classical sense, in that it does not rely on a diffusion process 
in velocity space induced by the discrete number of particles in the system. 


(C9) 


(CIO) 


C2 Linking the infall, drift and diffusion to the cosmic two-point correlations 

Up to this point we investigated the secular evolution of a given (phase averaged) halo, undergoing a given inflow and tidal 
field accretion history. Let us now invoke ergodicity so as to replace temporal averages by ensemble averages in Eqs. Ol-Ool. 
In doing so, we now try and describe a mean galactic halo embedded in the typical environment presenting the most likely 
correlations. This involves replacing ( )'[-with( )=F{ 1. Let us use Eg. ITot to expand Eg. <C10< . This yields: 


(D2(i,t))= x:k®kx: f r 

k n,n' ^ 


«'(0an(T))e*'"('' kdr) '*(!)!/;["'(I), 


(Cll) 


where an(t) = fln(f) + bn(t) corresponds to the coefficient of the total (self-consistent plus external) potential. If the first-order 
perturbations are stationary, let us write the two-point cross-correlation of the temporal fields, (an(t), (A) C[an, aj,/] (t — t) 

so that the integral in Eq. Eni may be carried as (assuming parity for the correlation function): 


note that when = s'^ = 0 this equation is conservative by construction. 
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POO 

/ C[an, an'] (AT)e'’‘"'^MAT = [k-w], (C12) 

Jo 

giving the temporal power spectrum evaluated at the temporal frequency, k • cu. Consequently, the diffusion coefficient becomes: 

(C13) 

k n,n' 

The same procedure may be applied to the other coefficient: 

(Di(I)) = Ek E k4 • cu]) , (C14) 

k n,n' 

while, for the secular correlation, Eq. <C8> : 

<DoW) = ^ /(S.)rdw + EkEk4 [^ ■ co]) , (CIS) 

where Pac" [re] is the mixed power spectrum given by (anCn) = ([aj + ^n]en) • Hence 

Pa"/[k ■ 07] = ((Ab + l)x(l) • (b* ® c) + Acx(l) • (c* ® c)) [k • 07]. (C16) 

Recall also that (given Eq. <19> and Eq. <491 1 


p^'"'[k • 07 ] = ((Af, + l)x(A^ + 1) ■ (b ® b*) + AaxA* • (c ® c*) + (Af, + l)xA* ■ (b ® c*) + Aax(A^ + 1) • (c ® b*}) [k • 07]. (C17) 

where Aj, and Ac involve K and therefore the secular distribution function, F, via Eq. <161 . Recall that Af, and Ac involve 
(1 - K)^^, which reflects the fact that the perturbation is dressed by the self-gravity of the halo. Eq. <C7I . toeether with Eos. <C13I - 
<C14< and Eq. <C17I provides a consistent framework in which to evolve secularly the mean distribution of a galachc halo within 
its cosmic environment. Note that it is possible via Eq. <511 to apply non-linear correchons to the induced correlation within R2oo- 


APPENDIX D: PERTURBATION THEORY TO HIGHER ORDER 
D1 Perturbative dynamical equations 

In this sechon, we 'solve' the dynamical equation to order n, which will allow us in the next sechon to present the N-point 
correlation to order n. 


ik 


Dl.l Perturbation theory to all orders 

Recall that for n ^ 2, ff''^ (I, t)obeys Eq. <34< . Given Eq. <15< it follows that 

= E / dTexp(7k- a7[T- t]) V) + <5? I’qM] f dli/7j^‘''(l)i/7jf'*(l) ^ ■ 

E E / dTexp(ik-a7[T-t])[fl®(T)+4fcq(T)] ({2n)^ ( dl|i/7[‘ll(w,I),/("^'')(w,I,t)| i/7[Pi*(I)'j , (Dl) 

k=l q,k \ J / 

where the first term in Eq. <D1< corresponds to the usual self-gravity coupling at order n, and the sum corresponds to the feed 
of lower order potential coupling into the n* order equation. Here (I, t) obeys 

/k’'^(I/f) = E /dTexp(zk-a7[T-f]) X ^ ■ty^ip^^\a^f\T)+S^bci{T)]+ 

q •' L 01 jt=i 

E [ dTexp(;k-a7[T-f])cq(T)4£7k''‘’'. (D2) 

q J 

Note that the response in Eq. <D2< is, as expected, out of phase with respect to the potential excitation, a-a{T) because of inertia 
(hence the modulation in exp ( 7 k ■ w{t — t))). Now, ton^^' order Eq. <DH . <D2< maybe rewritten formally as (using the contraction 
rule Eq. 


-f 


= Ki • a'") + K2 • ( E + 4^] ® + 4^]) + • • • +K; ■ ( 

V<i+' 2 =n / y, 

® c) + ■ • • + Qy ( E + ill’ll ® 


E (8)[a‘''^+4b] •■•+Kn-(g)]aW+b]-f 

'U-Hi,=17 j j n 


Q 2 ■ a' 


C + • • • + Qn 


aW +b] 


n-1 


) C, 


(D3) 







30 Christophe Pichon and Dominique Aubert 





Figure Dl. diagrammatic representation of the expansion to third (top diagram, corresponding to Eq. IDSl i and fourth (bottom diagram, given 
Eq. ID6I ') order; again (see Fig. |5| for details) the coupling of the tidal interaction (through the bn coefficients) and the incoming infall (expanded 
over the Cn coefficients) yields the coefficient of the response inside R200- The coupling is achieved via the operator K; and Q; as explained in 
Fig. 1021 : the curly brace in front of the diagrams account for the number of such diagrams entering the expansion, corresponding to the 
permutation of the input (recalling that the order matters). Note also for each branch the sum of the order of the sub branch correspond to the 
order of the expansion. 


where the kernels Kj and K 2 are given by Eqs. <41> -<42>. while Kn,n ^ 2 obey formally: 


(K, 




E 

^^2n-l + 


,„[Ti-f,T2-Ti,---,Tn-Tn_i] = [27r]3^ / dIexp(;k-a;[Ti-f]) |[exp(tki • a7[T2 - Ti]) • 

V T ki+k2=k 


exp((k2n-3 ■ a3[rn - 




ri"'’ - 


(D4) 


Note that the n* order Kernel involves "only" one integral over action space, but n coupling in configuration space and n + 1 
time ordered instants (f, Ti, ■ • • Tn). Note also that Eq. ID1I implies that secular perturbation theory accounts for both the rate of 
change in frequency of the system, via 3” cu / the rate of change in equilibrium via / 3F but also the rate of change in the 
incoming flow via 3"c7'[Pl''’/3I". Note finally that the relative phases (causality) are accounted for via the ordered time integrals. 
For instance Eq. <D3t reads to third order: 

= Ki • + K2 • + b] ® ® +b] j+ 

K 3 • +b] ® [a^^^ +b] ® [a^^^ +b]^ + Q 3 • ([af ^ +b] ® [af^ +b] ® + Q 2 ■ a^^^ ® c, (D5) 

and is illustrated in Fig. <DH together with a^^^: 

a(^) = Kj • a^^) + K 2 • ^a^^^ ® [a^^^ + b] + [a^^^ + b] ® a^^^ + a^^^ ® a^^^ j + 

K 3 • ^a^^) ® [a^^^ + b] ® [a^^^ + b] + + b] ® ® [a^^^ + b] + [a^^^ + b] ® [a^^^ + b] ® a^^^^ + 

K 4 • [a^^) + b] ® [a^^) + b] ® [a^^^ + b] ® [a^^^ + b] + 

Q 2 • ^a^^^ ® cj + Q 3 • ® [a^^^ + b] ® c + [a^^^ + b] ® a^^^ ® cj + Q 4 • [a^^^ + b] ® [a^^^ + b] ® [a^^^ + b] ® c. (D 6 ) 

Note that Eq. <D6> depends recursively on Eq. <D.5> and both depends recursively on Eq. <441 and i43l . When the recursion is 
carried through, (see Fig. <D2» the expected relative complexity of the non-linear evolution appears clearly. 
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D1.2 Reordering to higher order 

In the main text, we give in Eq. 1531 and above the first- and second-order reshuffling of the perturbahon in b and c. Similarly, 
the third-order term reads in terms of products of b and c as 

• b ® b (g) b -b Accc • c ® c ® c -f • b ® b ® c + A^^i, • c ® c ® b -b 
• b ® c ® b + • c ® b ® b -b Af,^,. • b ® c ® c -b A^f,^ • c ® b ® c, 

where (following the same convention as in the main text for the brackets) 

Ahbb = o Ki' -b o [K'/, o K'/] + o [K^ o K'/, K'/], 

O Qi + o { [Q;, o Q; + o [Q;, 1]] + [K^ o + Q' o [Q[, 1], Q^]} -b o [Q[, Q[,l] + 
Q^o[K^oQ;+Q^oQ;,1],1], 

Abbe = ^3 ° [Ki, K'{, Q'l] + o [K'{, K'{, 1] + o [K'{, o [K'{, Qi] + o [K'{, 1]], 

Abcb = K^o[Ki',Qi,Kl']-bK^o[K^o[Ki',Qi]+Q^o[Ki',ll,Kl'], 

Accb = ^3 ° [Qi. Qi. Ki'] + o [Q'l, K'2 o [Q[, K"]] , 

Acbb = K^o[Q",K'/,K"]+K^o[q;,K^oK'/], 

Abec = K^ o [Kf Q'l, Qi] + o [Kf Q[,l] + o [Kf o -b o [Q^, 1]] + o [K", o -b o [Q^, Ij] + 

Q^o[K^o[K" Qi]+Q^o[K'/,l],l], 

Aebc = K^ o [q;, K'/, Q[] + o [Q", K'l, 1] + o [K^ o [Q^, K"], K'/j -b o [Q^, o [Kf Qi] + o [K", 1]] -b 

Q^o[K^o[q;,K"],1]. (D7) 

Generically, after reordering, Eq. <53> becomes 


ffft) = 


(il ® ■ 


it), 


l ii,---inG[b,c] 


ft fT„-l 

£ / dTi---/ dTn X] [Ai ■i„]p,qi...q„(f-T'l.'"Tn-Tn_i)[ii]q^(Ti)...[in]q„(Tn), 

il,-■i„G[b,c]'^-°° qi- qn 


which involves 2" terms. Here [Aj^. ..j^] 


(D8) 


p,qi ■ q,. 


( 01 , • • • 6n) is some linear tensor of order n -b 1 which returns the n*^'order response 


to the excitation b;(6),cy(6) at various times 6i,02 • • • 0p- Note that it involve the equilibrium distribution function, Eg and its 
derivatives with respect to the actions, I, together with the properties of the basis function. 


D2 The N-point correlation function 

In the main text, we presented the calculation of the two-point correlation of the fields within the R 200 sphere. More generally 
we are interested in the N-point correlation of, say, the density (at various times): 

00 

Cn = (p(h)| 0 ( 3 : 2 )---| 0 (xn)) = E E 

n=N Pi+P 2 + ---PN=n 

n=N pi+P2+---pN=n qi,--',qN 

Now, solutions to the n'^' order perturbation theory are given by Eq. <D8> . It follows that 


(D9) 


£ [dfteU...:] (Ti,ei,i,---, 01 , 

L ^lJqi,q,i...q,,^ 

f d^’^e (hM.eypi • • - .^pn.pn) X 

J L PNJq„,q J...q 


(«q'['^(Tl)'"flq^[r'(TN)} = E ■■■ E 

il/- ■ -ipi ^ [b,c] il,■ ■ -ipj^ e [b,c] qi,i'''4 pn'Pn 


([il]qu(®l,l) ■ ■ ■ [ipilqi.pi (®PN,l) ■ ■ ■ [ipNlqi.pfj (01 ,Pn) ■ ■ ■ [ipN]qpNPN(®PN,PN)) • 

If the perturbation is a centered Gaussian random field, Wick's theorem states that: 

([ii]qi,i(®i,i) ■ ■ ■ Ppilqi.pi (®pn,i) ■ ■ ■ [ipNlqi.pfj (^i,pn) ''' [ipNlqpijppj (®pn,pn)) = 

E n( [il]qu(^l,l)[ipi]qi,Pi (^PN,l)) ■ ■ ■ ([ipN]qi,pp,(^l,pN)PpN]qpMPN[®PN'PN)) ' 


(DIO) 


(Dll) 


all permutations 

Putting Eqs. IDlOt - lDllI into Eq. <D9t yields formally the N-point correlation function to arbitrary order. A special case in given 
in the main text corresponding to third order expansion of the two-point correlation, Eq. 1.541 . The N-point correlation of other 
(possibly mixed) moments of the distribution fimction may be computed following the same route. 
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Figure D2. Reordered diagrammatic representation to second-order (first-order included) of the expansion given in Eqs. 1431 - 1441 . This time only 
fcn and Cn are inputs. The closed loop accoimts for the self-gravity and represents (1 — The thin loop traces the fact that the perturbed 

potential contributes also directly to the second-order term via ai + b (see Eq. 1441 for details). Note that in each diagram, each oblique line 
represents a sum over k and a time integral. The dashed line stands for infall coupling, while the thick line stands for the tidal coupling. 


D2.1 Synthetic hierarchy 

D3 Perturbation theory in the complex Fourier plane 

Let us close this appendix by a presentation of the perturbative solutions in the complex Fourier plane. In frequency space, 
Eq. reads: 


3 ( 2 ) 


• a; — a; 


(") = E If • 

[2nfY^ / dl E / ^ 

k J qi,q2 J K • 07 - W 

If • + Sq2("')] +tr^'|‘*''(I)Cq2[07'] 


E 

kq+k2—k 


T 1^2 


kj • 07 — (o; — cij') 

Following Eq. <451 . let us also define in frequency space the contraction rule: 
(K-Z)p(o7)4£4,q(o;)Zq(o;), 

q 

(note that Eq. <D13> only involve a sum and no integral) and the higher order contraction rule (cf. Eq. <460: 

(Kn-Zi®---®Z")p(a;)= [ dw 

qi, • q^-f 

The operator, Kn[oti, 022 , • • •, Wn], obeys 




Kn ■ Z^ ® ■ ■ ■ ® Z^)p (07) — E f dcVi ' ' ' [ dt77nKp,qn .q„ (U7i, • • •, OJn)‘^o (o7 — E (^l) ' ' ' 2q^ (^n) ■ 

qi, ■q„-' J ,=1 


(Kn)p,qi,q2 "vijn [^1/^2/■ ■ ■ / Wn] — [27r]^ E / dl 

k •’ 


E 

ki+k2=k 


k;[ • 07 — 072 


k 2 «—lTk 2 n—kn 


k • 07 — Otj 

dF 


k2n-3 • 07 — 07n 91 






.♦IT' 


♦I". 


(D12) 


(D13) 


(D14) 


(D15) 


























APPENDIX E: OTHER COSMOLOGICAL PROBES 

In this appendix, we discuss other non-linear statistical 
probes of the cosmic environment of haloes, expanding over 
Section lSrl 


El Dark matter disintegration 

It has been claimed that dark matter could be made of neu- 
tralinos which can be traced indirectly via their disintegra¬ 
tion signature, which scales l ike th e square of fhe local dark 
matter density iStnehr et alJ i20n,'^i i. The total number of 7 
photons received during integration time, reads 

Na(n,t,) = 

X 

4 J dn J drpl^ir), (El) 

where Dgff is the effective size of the telescope, AQ the an¬ 
gular resolution of fhe telescope, and NcontiEj) is the num¬ 
ber of continuum photons and (crv) is the continuum cross- 
section of neutralinos of mass, The integral accounts for 
measured flux of 7 photons arisi ng from neutralino s desin- 
tegrating in the direction, O. tSee lStoehr et alJ i2n03h for de¬ 
tails about the computation of Ncont, Deff Since the 

Na(n, T) scales like the line integral of the square of the den¬ 
sity along the line of sight, it is straightforward to propagate 
the statistical properties of the density fluctuations to that of 
Na. For insfance, the cosmic mean will scale like 

(Nannih(n)) = W;, J dr{puM{r)f+ Wx J dr{Spl^{T)) ,(E 2 ) 

where 

(^Pdm(J')) = (E3) 

n,n' 

Hence we expect an excess of annihilation because of fhe po¬ 
larized clumps within the halo. Similarly, we may predict 
the angular correlation function, or the related variance as 
a function of smoothed angular scale as 

(JNa(n)<5Na(n')) =</ drdr'{SplM{n,r)SplM(n'y)) 

= E // drdr' X 

ni,n2,n3,n4 

{r, {r', o'), (E4) 

where ^Wannih{^) ^annih(^) ~ (^annih{^))- Note that we 
assumed here that the resolution of fhe telescope was effec¬ 
tively infinite (i.e. AQ —> 0 in Eq. lEIH . Now we may rely 
on Wick theorem to express the four-point correlation en¬ 
tering Eq. IE4> as products of two-point correlations. Calling 
6a = a — (a), we have {Scin^San 2 ^an 3 ) = 0 , and 

If the infall is statistically isotropic, Eq. <E4> may be averaged 
over the direction, Q and reads: 

{SKnnih{Ci)SN^rmih{Cl')}n = ■ fi'], (E 6 ) 

I 
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^^arinih _ 

~ 

i 

Note that the geometric factor, only depends on the 

basis function, pi"! (r) and possibly the resolution of fhe tele¬ 
scope if it is not assumed to be infinite: 

^Lh= E /dnyf*(n')/dn'yf*(n')x 

JJ drdr'p^’^^\r,(n)p^’^^\r,a)p^''^f(n',r')p^’^^fn',r'), (E 8 ) 

given that p["l(r) = u'gyr)Y'^{Cl) and given the 

properties of spherical harmonics, the integral 

j dClY’f^{n)Y”f{n)Yff{n')Yff{Cl')dCl can be re ex¬ 
pressed iteratively in terms of Clebsch-Jordan coefficients. 
Ensemble average and comparison with the observation is 
possible at the high £ limit corresponding to the small-scale 
structure of the dark matter halo, for which we may expect 
independent angular regions of the Galactic halo to be 
representative of an ensemble average. 


E2 Bremsstrahlung X-ray emission of stacked haloes 


Assuming that the gas traces the dark matter, we may re¬ 
produce the thought experiment of Section lEll though the 
ensemble average is constructed while staking projections of 
haloes on the sky rather than in a galactocentric framework. 

The emissivity per imit vol ume at frequen cy v, £v{t), for 
a hydrogen plasma is given by iPeacoc 3 EH) 


£„(r)drdi/ = 


ex»e(r) 

V^) 


1 + logio 


HTeit) 


hv 


exp 


hv 


ksTeir) 


drdv. 


(E9) 


where Tg is the temperature in Kelvin, v the frequency in 
Hz, kg the Boltzmann constant, h the Planck constant, and 
ex = 6.810^^^ for an emissivity in Wm^^Hz^^. Let us assume 
here that the cluster is isothermal, hence the variation of Tg 
with z are neglected compared to that of Ug squared.Let us 
also assume that L/M is the mass to light ratio of the cluster 
is constant. Hence the emissivity per unit surface, ay, is given 
by 


(7-4R) = j dz{L/M)^ey{R,z)^Wx JdzplyR,z). (ElO) 


Hence, taking an ensemble average yields: 

(c7,(R)) = WxJ dz(pDM)"(R,z)+Wx J dz{Sply{R,z),(mi) 


where 


(Pdm)(R/Z)= (E12) 

n,n' 


The two-point correlation of the cosmic fluctuation of 
the emissivity is given by 


where 


13 


this is a better approximation than for the SZ effect 
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{SaPR)5a4R')) _ JJ {Spl^{K,z)5pl^(R',z')) 

Av{R)P Av(R)P 

where 

(^|0dm(R.z)‘5Pdm(R^z')} = E (^nian^an^and X 

ni,n2,n3,n4 

p[''il(R,z)pl"J(R,z)p[''3l(R'^z')p[n4l(R^^2'). (E13) 

Note the cancellation of the dependence on Wx (hence T or 
M/L) in Eq. <E13> . Relying again on Wick's theorem, Eq. <E5> . 
we may express the four-point correlations as a products of 
known (cf. Eq. <E5I 1 two-point correlations. 


E3 Galactic halo's ellipticity 

More generally, let us consider a problem which depends 
non trivially on the perturbed distribution function, e.g. the 
ellipticity, Ch/ of fho departure from sphericity of the sub¬ 
structures induced by the environment around a given halo. 
The ellipticity is defined as 

ch = - 1 = S{Sp{i )), with {A,} = Eigenval(lH) ,(E14) 

L^i ^7 

and 

%!,;=/ diSp{t)xiXj / / drpNFw(i-), (E15) 

^^200 / -7 ^^200 

(so that Aj is the largest eigenvalue of Ih and ch = 0 if the 
halo's perturbation is spherical). Since we know the statisti¬ 
cal properties of Sp{t), we may predict the statistical proper¬ 
ties of ch- In practice, assuming 5 is a well-behaved function 
of its arguments, we may Taylor-expand Ch with respect to 
Sp as : 

= E ■ • • l^PPn) - (<ip(r„))] (E16) 

Note that the derivative in Eq. <E16t is a Erechet functional 
derivative, so that the dot involves an integration over r. 
Hence the ensemble average, (ch) will involve N-point cor¬ 
relations, and reads 

<®h)=EE ■ P^'^HP ■ ■ ■ P^'"\in), (E17) 


where, once again, we may rely on Wick's theorem to reex¬ 
press (fl,j • • • a A as products of two-point correlations. Since 
the relationship between the density perturbation and the el¬ 
lipticity is not linear, we expect a non zero ellipticity on aver¬ 
age. 

Note that in principle, we may reconstruct the full 
PDE of e. Formally, calling z= (e, fl 2 /• • ■ ^n) (so that z = 
( 6 (^ 1 /• • ■ /fln),fl 2 • • ■ / fln) = g{ai, ■ ■ • ,fln)), inverting for z as 
a function of {«, } (provided the ellipticity is not degenerate 
in flj), and marginalizing over the other coefficients yields: 


PDF(e) = Jda2 - dfl„PDF (g-i(z)) / 


dz 


Now in practice, Eq. <E17I might not be the simplest pro¬ 
cedure to compute (cr)/ and monte carlo resimulahon may 
turn out to be more practical. 


E4 Metal lines in QSO DLA systems 


Let us finally consider a more convolved observable, which 
will depend on both the clump distribution within the 
haloes, but also on their velocities. 

In the red part of a high resolution spectrum of quasars, 
groups of absorption features are found, corresponding to 
the physical situation where the light emitted by the quasar 
is partially absorbed by the metal-rich^^ clumps which the 
line of sight happens to intercept. Formally, the normalized 
flux in a QSO is proportional to minus the log of the optical 
depth along the line of sight. Th e optical depth in the metal 
transition is jPichon et aljoooilh : 


t{w, R) 


cctq /•+” nz{v,R) 
H(z)\Ai J-oo b{v,R) 


exp 


(w — V — Vz(v, R))^ 
b(v, R )2 


dv, 


(E18) 


where c is the velocity of light, uq is the metal absorp¬ 
tion cross-section, Tf(z) is the Hubble constant at redshift z, 
nz{v,R) the ionized metal number density field, b{v,R) the 
Doppler parameter (accounting for the thermal broadening 
of the line), and Op (o, R) is the peculiar velocity, at impact pa¬ 
rameter, R from the centre of the cluster. The observed nor¬ 
malized flux, F, is simply F = exp(— t). If we assume here 
again constant biasing, so that «2 Pdm- This assumption 
may be lifted once the identification of virialized substruc¬ 
ture described in Section l5.3.1l is carried through. The two- 
point correlation of the optical depth fluctuation will involve 
statistical properties of both the density and the velocity field 
in a non trivial manner. 


1 

{Tp{w,R) 


R)St{w', R)) , 


(E19) 


with 


(iT(H;, R) = T(zi', R) — {t){w,R) . (E20) 

Note that the distance to the halo center, 
R= b(cos[!9{,],sin[!9f,]) still occurs in Eq. IE20> . Since we 
do not know in general the impact parameter of the line of 
sight with respect to the halo center, let us marginalize over 
its a priori probability distribution, which we may infer from 
e.g. the PT model (which at these scales corresponds essen¬ 
tially to the autocorrelation of the unperturbed universal 
halo profile). Given that we consider systems at the redshift 
of a damped Lyman-a, we may assume that we fall close to 
a galactic structure. Calling pi,{b,z,M)dbdz the probability 
of a given point in space to be at a distance, b within db of 
an object of mass larger than M, which is at redshift z within 
dz, we may construct the weighted sum : 

^00 nCO n27T 

Ct(Ah;) = db dz d&i,pi,{b,z, M) x (E21) 

Jo Jo Jo 

{St{w, b cosji?;,], b sin[!9j,])i5T(ro -i Aw, b cos[&i,],b sin[!?f,] ))jj, . 

This quantity may now be compared to the observable. Let 
us assume some equation of state for the metal phase, so that 


Since we make predictions at lower redshift we need to concen¬ 
trate on metals such as Mgjj, or Fen which are found typically at 
redshift z ^ 1.5 in the visible 
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b{R,z) = bg (p(R,z)/ py. Eq. IE 181 may then be written for¬ 
mally as bT(w,R) = Tlbp(v,R),Vz(v,R)]. Let us Taylor ex¬ 
pand this expression in the neighborhood of the mean den¬ 
sity fluctuation as: 

[(ip(ri) - (^p(ri))] • • • Ibvz(rn) - (^Uz(r«))] • (E22) 

Again the derivative in Eq. IE22I is a functional derivative 
(cf. Section lE3l . Eqs. IE21I - IE22I together with Eq. 1291 yield 
the expected correlation as a fimction of fhe statistical envi¬ 
ronment. 



